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The  injection  of  freshwater  into  a  saline  aquifer  can,  under 
favorable  conditions,  result  in  displacement  of  the  saltwater  and  crea- 
tion of  a  distinct  body  of  freshwater  called  a  freshwater  lens.  It  may 
be  desirable  to  do  this  either  as  a  method  of  storing  the  freshwater 
for  future  use  or  in  order  to  form  a  barrier  against  saltwater  intrusion. 
In  either  case,  it  is  important  to  know  beforehand  whether  the  lens  will 
be  stable  or  the  injected  water  will  be  lost.  The  purpose  of  this 
study  is  to  develop  methods  for  analyzing  the  hydrodynamic  factors 
involved  in  the  formation  and  maintenance  of  such  lenses. 

A  series  of  laboratory  tests  using  vertical  Hele-Shaw  models  was 
performed  to  determine  the  shape  and  rate  of  movement  of  the  initially 
vertical  interface  between  two  fluids  of  different  density.  The 
results  of  these  tests  were  compared  with  the  predictions  of  two  formulas 
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proposed  by  previous  investigators.  It  was  found  that  for  relatively 
thick  aquifers  it  is  not  correct  to  neglect  the  vertical  components  of 
flow  as  has  formerly  been  done. 

An  analytical  method  is  presented  for  the  location  of  the  boundary 
of  a  steady  freshwater  lens  in  a  coastal  aquifer.  A  vertically 
integrated  potential  function  is  used  to  superimpose  the  flow  from  the 
injected  well  on  the  pre-existing  uniform  seaward  flow  in  the  aquifer. 
This  results  in  a  quasi-three  dimensional  solution  for  the  position  of 
the  coastal  interface. 

A  pair  of  coupled  ordinary  differential  equations  is  derived  which 
describe  the  steady  flow  of  freshwater  injected  into  a  leaky  saline 
aquifer  under  the  Dupuit  assumption.  The  non-linear  system  of  equations 
is  solved  numerically  by  isolating  the  unique  integral  curve  which 
results  in  the  only  physically  possible  boundary  conditions.  A  design 
chart  is  presented  by  which  the  maximum  possible  size  of  freshwater 
lens  can  be  determined  for  various  aquifer  conditions. 

A  three  dimensional  finite  element  model  is  developed  which  can  be 
used  to  locate  the  boundaries  of  a  steady  freshwater  lens  under  a  wide 
variety  of  conditions.  The  model  uses  isoparametric  elements  of  the 
quadrilateral  family  having  twenty  nodes  each.  This  permits  very  close 
approximation  to  the  curved  flow  boundaries  typical  of  freshwater  lenses. 
Inhomogeneous  and  anisotropic  aquifers  can  be  treated. 

A  finite  difference  model  developed  for  the  U.S.  Geological  Survey 
is  used  to  analyze  the  axisymmetric  flow  in  transient  freshwater  lenses. 
Simulation  of  lens  formation  in  a  variety  of  saline  aquifer  configura- 
tions indicates  that  previously  available  methods  of  analysis  which 
neglect  vertical  components  of  flow  may  be  reliable  in  very  thin 
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aquifers.  For  thicker  aquifers,  however,  vertical  flow  can  predominate 
resulting  in  a  much  more  rapid  rate  of  gravitational  segregation  of 
the  fluids. 
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CHAPTER  1 
INTRODUCTION 

Groundwater  flow  is  an  important  part  of  the  hydrologic  cycle  which 
nature  uses  to  even  out  the  temporal  and  spatial  variations  inherent  in 
rainfall.  It  has  long  been  man's  practice  to  take  water  from  groundwater 
reservoirs  with  the  assumption  that  the  water  removed  will  be  replaced  by 
natural  processes.  In  many  cases,  this  arrangement  works  well,  with 
nature  putting  water  into  the  ground  and  man  taking  it  out.  In  some 
areas,  though,  the  increased  use  of  groundwater  supplies  coupled  with 
decreased  natural  recharge  due  to  man-made  changes  in  the  land  surface 
is  either  draining  the  freshwater  aquifers  or  opening  them  to  saltwater 
intrusion.  In  recognition  of  this,  efforts  have  been  made  to  augment 
natural  recharge  by  such  means  as  spray  irrigation,  water  spreading, 
seepage  ponds  and  deep  well  injection.  This  signals  a  new  conception 
of  aquifers  as  a  form  of  underground  storage  which  can  be  used  in  much 
the  same  way  as  surface  reservoirs.  As  an  extension  of  this,  it  has  been 
proposed  to  make  use  of  the  storage  capacity  available  in  saline  aquifers 
by  injecting  freshwater  into  them  in  such  a  way  that  the  resident  salt- 
water is  displaced  and  a  lens  of  freshwater  is  created. 

The  variety  of  possible  flow  patterns  that  may  result  from  the 
injection  of  freshwater  into  a  saline  aquifer  can  be  viewed  as  a  contin- 
uous spectrum  bounded  by  three  idealized  extremes.  One  possibility  is 
that  the  freshwater  will  disperse  rapidly  and  lose  its  identity  as  a 
separate  flow,  merely  serving  to  dilute  the  resident  saltwater.  In 


another  case,  a  flow  pattern  may  develop  which  transports  the  freshwater 
away  from  the  point  of  injection  as  a  plume.  A  third  alternative  is 
that  the  freshwater  will  tend  to  remain  near  the  point  of  injection 
creating  an  identifiable  lens  of  freshwater  which  grows  or  diminishes 
accordingly  as  fluid  is  added  to  or  removed  from  it.  If  the  freshwater 
is  injected  in  anticipation  of  its  subsequent  retrieval  and  re-use, 
the  formation  and  maintenance  of  such  a  lens  adjacent  to  the  point  of 
injection  is  essential. 

The  major  factors  influencing  the  formation  of  freshwater  lenses 
are:  (1)  buoyancy,  (2)  mode  of  injection,  (3)  the  properties  of  the 
porous  medium,  (4)  dispersion,  (5)  flow  of  the  resident  saltwater. 
The  relative  intensities  of  these  individual  influences  determine  which 
of  the  previously  described  flow  patterns  will  predominate.  Favorable 
conditions  for  the  development  of  freshwater  lenses  require  that 
buoyancy,  dispersion,  and  preexisting  flow  be  relatively  weak  influences. 
In  previous  studies  of  freshwater  lenses  such  ideal  conditions  have 
been  assumed.  In  reality,  these  effects  are  never  entirely  absent. 
It  is  important,  then,  to  determine  at  what  point  their  effects  become 
substantial  and  to  have  techniques  for  including  their  influences  in 
the  analysis  of  freshwater  lenses. 


CHAPTER  2 
LITERATURE  SURVEY 

General 


The  analysis  of  freshwater  lenses  is  one  of  a  general  class  of 
problems  involving  the  stratified  flow  of  two  miscible  liquids  through 
porous  media.  It  is  closely  related  to  such  other  problems  of  that 
class  as  saltwater  intrusion,  interface  upconing  beneath  wells,  and 
deep  well  disposal  of  liquid  wastes.  A  wealth  of  published  information 
is  available  on  each  of  these  topics  but  much  of  it  has  no  direct 
application  to  the  study  of  freshwater  lenses.  Many  of  the  methods 
developed  for  locating  coastal  interfaces,  for  instance,  depend  on 
conformal  mapping  techniques  or  otherwise  on  the  theory  of  complex 
variables,  which  is  limited  to  two  dimensional  problems.  Likewise, 
much  of  the  work  done  with  deep  well  disposal  deals  with  the  movement 
and  dispersion  of  plumes,  emphasizing  the  rapid  departure  of  the  waste 
fluid  from  the  scene  of  injection.  Nonetheless,  some  of  the  work  in 
each  of  these  fields  can  be  useful  in  lens  analysis.  It  is  only  these 
studies  and  those  which  are  directly  concerned  with  fresh  or  hot  water 
lenses  that  will  be  mentioned  here. 

Studies  Dealing  with  Interface  Location 
The  rotational  movement  of  an  initially  vertical  interface  between 
two  fluids  of  different  density  was  studied  by  Gardner,  Downie  and 
Kendall  (1962)  because  of  its  applicability  to  recovery  processes  in 
petroleum  reservoirs.  The  initial  motion  of  the  interface  and  the 


long  term  motion  were  found  analytically  using  the  assumption  that  the 
porous  medium  was  thin  and  vertical  flow  could  be  neglected.  These  two 
solutions  were  then  joined  to  give  one  equation  which  was  assumed 
accurate  for  intermediate  times  as  well.  It  was  assumed  in  these 
derivations  that  the  interface  remained  straight  during  its  rotation 
toward  the  horizontal.  Laboratory  tests  using  models  of  square  and 
circular  cross-section  packed  with  glass  beads  were  performed  which 
seemed  to  agree  with  the  equation. 

Dagan  and  Bear  (1967)  treated  the  upward  movement  of  a  freshwater- 
saltwater  interface  beneath  a  partially  penetrating  well  using  a  pertur- 
bation expansion  to  linearize  the  interface  boundary  condition.  They 
checked  their  results  with  a  sand-packed  model  and  found  them  to  be 
reasonable  when  the  disturbance  of  the  interface  was  small.  The  range 
of  applicability  of  this  method,  however,  appears  to  be  quite  limited. 
Shamir  and  Dagan  (1971)  developed  a  transient  numerical  model  to 
predict  the  landward  movement  of  a  saltwater  wedge  in  a  coastal  aquifer. 
The  Dupuit  assumption  was  used  to  formulate  the  governing  equation  which 
was  then  solved  by  an  implicit  finite  difference  method.  Dispersion  was 
neglected.  The  numerical  results  were  checked  against  an  analytical 
solution  similar  to  that  of  Gardner,  Downie  and  Kendall  and  with  a 
Hele-Shaw  model.  The  comparison  with  the  analytical  solution  showed 
good  agreement  that  the  interface  moved  as  a  straight  line.  Comparison 
with  the  Hele-Shaw  model,  however,  was  not  quite  as  good.  The  authors 
attributed  this  to  the  fact  that  vertical  velocities  were  neglected  in 
both  the  analytical  and  numerical  solutions. 

Strack  (1976)  studied  the  influence  of  inland  production  wells  on 
the  steady  state  shape  of  the  saltwater  wedge  in  coastal  aquifers.  His 


analysis  makes  use  of  a  vertically  integrated  potential  function  which 
can  be  defined  so  as  to  be  single-valued,  throughout  the  aquifer. 
Lines  of  constant  potential  can  then  be  found  in  the  horizontal  plane 
which  lead  to  a  three-dimensional  description  of  the  interface  based  on 
the  Ghyben-Herzberg  principle. 

A  similar  problem  was  studied  by  Kishi  and  Fukuo  (1977).  They 
considered  the  effect  of  inland  pumping  on  the  saltwater  wedge  in  a  fan 
shaped  coastal  aquifer.  The  Laplace  equation  in  radial  coordinates  is 
integrated  vertically  to  make  it  two  dimensional.  Using  the  Ghyben- 
Herzberg  principle  the  dependent  variable  is  transformed  so  that  it  will 
be  single-valued  in  both  the  freshwater  region  and  the  interface  region. 
The  resulting  equation  is  solved  by  the  method  of  Green's  functions  to 
give  a  three-dimensional  description  of  the  interface. 

Studies  Dealing  with  Dispersion 

Gardner,  Downie  and  Wyllie  (1962)  considered  the  mixing  across  the 
interface  between  methane  gas  and  nitrogen  used  as  a  cushion  gas  in  a 
natural  gas  reservoir.  A  solution  to  the  convective  dispersion  equation 
in  radial  coordinates  was  combined  with  the  gravity  segregation  formula 
found  by  Gardner,  Downie  and  Kendall. 

Harleman  and  Rumer  (1962)  derived  an  equation  for  interface  rota- 
tion which  is  the  same  as  that  found  by  Gardner,  Downie  and  Kendall  for 
the  latter  stages  of  movement.  This  equation  also  is  based  on  the 
assumption  of  a  straight  interface.  The  same  report  also  deals  with 
longitudinal  and  lateral  dispersion  in  both  stationary  and  moving 
coastal  interfaces.  Laboratory  tests  using  sand-packed  models  to 
simulate  a  coastal  aquifer  showed  that  in  the  absence  of  tidal  fluctua- 
tions the  mixing  zone  was  quite  thin,  especially  near  the  toe  of  the 


interface.  An  equation  was  presented  which  predicted  the  thickness  of 
the  dispersion  zone  accurately  in  this  region.  When  changing  conditions 
in  either  the  saltwater  or  the  freshwater  caused  a  shift  in  the  equi- 
librium position  of  the  interface  the  effects  of  longitudinal  dispersion 
became  pronounced.  It  was  concluded  that  such  shifts  are  responsible 
for  most  of  the  dispersion  found  in  natural  coastal  aquifers. 

Wooding  (1963)  obtained  approximate  solutions  for  the  lateral 
dispersion  across  a  stationary  interface  by  transforming  the  governing 
equations  into  a  natural  coordinate  system  based  on  the  shape  of  the 
interface.  Using  order  of  magnitude  reasoning  these  equations  were 
reduced  to  the  form  of  Prandtl's  boundary  layer  equations.  For  these, 
the  usual  similarity  solution  was  obtained  which  relates  the  thickness 
of  the  mixing  layer  to  the  distance  from  the  stagnation  point  and  the 
slope  of  the  interface. 

Hoopes  and  Harleman  (1965  and  1967)  analyzed  the  dispersion  of 
tracers  injected  into  groundwater.  Their  work  consisted  of  analytical 
and  numerical  solutions  to  the  convective  dispersion  equation  and 
laboratory  experiments  with  a  sand-packed  model  to  verify  these  solu- 
tions. The  relative  effects  of  convection,  lateral  and  longitudinal 
dispersion,  and  molecular  diffusion  were  compared  for  the  case  of 
axisymmetric  radial  flow  from  an  injection  well  and  for  the  flow  between 
an  injection  well  and  production  well.  It  was  shown  that  the  effects 
of  lateral  dispersion  on  the  steady  flow  between  two  wells  was  almost 
negligible  except  near  the  line  joining  the  two  wells. 

Studies  Concerned  Directly  with  Freshwater  Lenses 
Glazunov  (1967)  proposed  to  use  a  cluster  of  three  to  five  wells 
operated  in  various  patterns  of  injection  and  production  in  order  to 


form  freshwater  lenses  of  the  desired  shape.  His  analysis  neglects 
both  dispersion  and  gravity  segregation,  treating  the  convection  of  an 
interface  in  two  horizontal  directions  by  a  method  introduced  by  Muskat 
(1946).  In  this  method  the  interface  is  represented  as  a  material 
surface  which  is  transported  along  the  streamlines  by  the  flow  while 
exerting  no  effect  on  it  as  a  boundary. 

Esmail  and  Kimbler  (1967)  studied  the  effects  of  radial  flow 
geometry  on  gravity  segregation  and  longitudinal  dispersion  in  labora- 
tory models  composed  of  glued  sand  grains.  The  solution  to  the  convec- 
tive  dispersion  equation  given  by  Gardner,  Downie  and  Wyllie  was  verified 
for  cyclic  storage  in  freshwater  lenses  for  several  cycles  of  injection 
and  production  in  the  model.  A  dimensionless  parameter  analysis  of  the 
influence  of  radial  geometry  and  interfacial  mixing  on  the  rate  of 
gravity  segregation  was  conducted  resulting  in  an  empirical  equation 
for  the  rate  of  interface  rotation.  It  was  found  that  there  is  inter- 
action between  the  processes  of  dispersion  and  interface  rotation  which 
retards  gravity  segregation.  The  authors  concluded  that  the  horizontal 
projection  of  the  interface  is  proportional  to  the  square  root  of  the 
density  gradient  across  it  and  that  gravity  segregation  in  their  model 
had  more  effect  on  the  recovery  efficiency  than  did  dispersion.  Under 
these  conditions,  interface  dispersion  may  actually  improve  the  recovery 
efficiency  by  retarding  gravity  segregation. 

Kumar  and  Kimbler  (1970)  continued  the  research  with  glued  sand 
models,  finding  that  the  recovery  efficiency  improves  after  several 
cycles  of  injection  and  recovery. 

Whitehead  (1974)  developed  a  computational  method  for  predicting 
the  recovery  efficiency  of  a  single  well  cyclic  storage  system  in  thin 


horizontal  saline  aquifers.  His  method  was  based  on  the  empirical 
equation  of  Esmail  and  Kimbler.  Whitehead  verified  the  results  of 
his  computational  scheme  using  the  glued  sand  grain  model  and  found 
them  to  be  accurate  to  within  10%  for  three  storage  cycles.  The 
recovery  efficiencies  achieved  in  the  laboratory  tests  were  in  the 
range  of  70-100%.  It  was  concluded  that,  given  favorable  subsurface 
conditions,  the  use  of  saline  aquifers  for  freshwater  storage  could  be 
economically  feasible. 

In  a  series  of  reports  {Agrawal ,  1975;  D'Amico,  1975;  Tate,  1976; 
Kimbler  and  Whitehead,  1977)  the  results  of  Whitehead  were  extended  to 
include  the  effects  of  aquifer  dip  and  viscosity  difference  between  the 
resident  and  injected  fluids.  These  studies  also  used  glued  sand 
models.  The  conclusions  reached  from  these  experiments  were  that  the 
best  conditions  for  cyclic  storage  occur  when  the  saline  aquifer  is 
horizontal  and  the  injected  fluid  has  the  same  viscosity  as  the  resident 
fluid. 

It  is  widely  recognized  that  a  pre-existing  uniform  flow  in  the 
saline  aquifer  is  detrimental  to  the  recovery  efficiency  of  a  lens 
storage  system.  Langhtee  (1974)  proposed  to  surround  the  storage  area 
with  bounding  wells  which  would  be  pumped  in  such  a  way  that  a  stag- 
nation zone  would  result  inside  the  boundary.  The  required  pumping 
rates  for  these  wells  were  found  by  superimposing  the  streamline  and 
equipotential  pattern  for  the  individual  wells  and  solving  for  the 
pumping  rates  which  would  give  the  best  least-squares  fit  to  the  desired 
piezometric  head  on  the  boundary. 

Molz  and  Bell  (1977)  also  considered  the  use  of  boundary  wells  for 
head  gradient  control.  They  used  a  linear  programming  technique  to  find 


the  optimum  set  of  bounding  wells  which  would  give  the  desired  stagna- 
tion in  the  storage  area  with  the  least  expenditure  of  energy  for 
pumping. 

The  shape  and  movement  of  naturally  occurring  freshwater  lenses  on 
Grand  Cayman  Island  were  studied  by  Childley  and  Lloyd  (1977).  They 
used  a  modified  Ghyben-Herzberg  ratio  based  on  field  measurements  of 
salinity  distribution  in  the  mixing  zone  of  the  lens.  The  long  term 
transient  effects  of  seasonal  variations  in  recharge  rate  were  modeled 
numerically  using  the  idea  of  a  succession  of  steady  states.  Comparison 
of  the  numerical  results  with  records  of  actual  groundwater  conditions 
on  the  island  indicated  that  their  model  was  useful  in  predicting  long- 
term  changes  due  to  water  supply  wells  and  changes  in  recharge  condi- 
tions. They  concluded  that,  while  it  may  take  many  years  for  a  pumped 
lens  to  reach  equilibrium,  an  overpumped  lens  can  collapse  locally  in 
a  matter  of  a  few  months. 

Smith  and  Hanor  (1975)  report  on  a  field  test  of  the  freshwater 
lens  storage  concept.  A  10  centimeter  fully  penetrating  well  was 
drilled  into  a  confined  saline  aquifer  36  meters  thick  to  be  used  for 
injection  and  retrieval  of  freshwater.  After  the  construction  of  this 
well  it  was  discovered  that  a  uniform  flow  of  approximately  15  centi- 
meters per  day  existed  in  the  aquifer.  Recognizing  the  unfavorable 
conditions,  the  investigators  proceeded  with  two  tests,  each  using  about 
757,000  liters  of  freshwater.   In  the  first  test,  the  water  was  stored 
for  two  hours  and  41%  was  recovered.  In  the  second  test  the  storage 
period  was  six  days  and  the  recovery  efficiency  was  24,5%. 

Werner  and  Kley  (1977)  conducted  a  field  test  of  the  storage  of 
heated  water  in  a  water  table  aquifer  3.35  meters  thick.  A  total  of 
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430  cubic  meters  of  water  at  a  temperature  of  45°C  was  injected  and  the 
movement  of  the  heat  was  monitored  with  buried  temperature  sensors  for 
a  period  of  64  days.  The  heated  water  remained  close  to  the  injection 
wen  throughout  the  test  period  but  no  attempt  was  made  to  recover  it. 

Molz,  Warman  and  Jones  (1978)  reported  a  field  test  with  a 
partially  penetrating  well  in  a  confined  aquifer  21  meters  thick.  A 
group  of  14  observation  wells  were  constructed  to  monitor  the  movement 
of  the  hot  water  injected.  In  this  test,  storage  of  7570  cubic  meters 
of  hot  water  for  36  days  resulted  in  an  energy  recovery  factor  of  0.69 
which  was  considered  promising. 

Papadopulos  and  Larson  (1978)  applied  a  finite  difference  numerical 
model  developed  for  the  U.S.  Geological  Survey  (INTERCOMP,  1976)  to 
simulate  the  test  results  of  Molz,  Warman  and  Jones.  The  predicted 
recovery  efficiency  was  75%  versus  a  value  of  69%  measured  in  the  field. 
Considering  the  accuracy  of  the  measured  aquifer  properties,  this 
prediction  was  considered  quite  good. 


CHAPTER  3 
BASIC   EQUATIONS  AND  FLUID  PROPERTIES 

Basic  Equations 
Darcy's  Law 

The  groundwater  flow  to  be  dealt  with  in  this  report  is  limited 
to  laminar  flow  in  porous  media  with  only  primary  porosity.     The 
equation  of  motion  for  this  flow  is   Darcy's  law  which  can  be  expressed, 
in  the  cartesian  coordinate  system,  as 


\  -  ^xx  ^x  '  hy  'y  '  hz  \  (3.1a) 


^  ~-  V  Jx  '  ^yy  Jy  '  V  ^z  (3. lb) 

^z  -  hx  \  ^  hy  Jy  '  hi  ^  (3-lc) 

where 

q^'  Py.  q^     =     components  of  the  discharge  velocity  in  the 
x,  y,  and  z  directions,   respectively 
K^j     =     components  of  the  hydraulic  conductivity  tensor 
for  (i,  j     =     x,  y,   z) 

J       =     _  M 
x  9x 

•      J      =    .M 
y  3y 
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<})      =     piezometric  head 
If  it  is  assumed  that  the  coordinate  axes  are  parallel   to  the  principal 
directions  of  hydraulic  conductivity  then  Darcy's  law  can  be  simplified 
to 


^  XX  dx    '        "^yy  9y  'J   "^zz  3z  ''^•'^' 


/^     /\ 


where  i,  j,  and  k  are  unit  vectors  in  the  x,  y,  and  z  directions, 
respectively.  This  assumption  is  not  a  major  restriction  to  the  use  of 
equation  (3.2)  because  in  practice  it  is  usually  difficult  to  distin- 
guish any  directional  hydraulic  conductivities  except  the  vertical  and 
the  horizontal.  If  the  soil  is  isotropic  or  only  one  dimensional  flow 
is  being  considered  equation  (3.2)  reduces  to 

q  =-K  V<t>  (3.3) 

The  value  of  the  hydraulic  conductivity,  K,  or  of  its  directional 
components  depends  on  the  properties  of  both  the  soil  and  the  fluid. 
Since  the  properties  of  saltwater  and  freshwater  are  somewhat  different 
it  may  be  necessary  in  the  general  case  to  consider  the  hydraulic 
conductivity  as  a  function  of  the  salt  concentration.  The  flow  resis- 
tance of  the  soil  is  then  described  in  terms  of  its  intrinsic  permea- 
bility, k,  which  is  related  to  hydraulic  conductivity  by 

K  =  kg/v  (3.4) 

where  v  is  the  kinematic  viscosity  of  the  fluid  and  g  is  the  accelera- 
tion due  to  gravity. 
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The  piezometn'c  head,  <t>,   appearing  in  equations  (3.2)  and  (3.3)  is 
defined  as 


(t>     =     z  +  P/y 


(3.5) 


which  depends  on  the  unit  weight,  y,   of  the  fluid  concerned.  Since  y   is 
also  a  function  of  salt  concentration  it  is  often  more  convenient  to 
deal  directly  with  pressure,  P,  in  the  equation  of  motion. 

Continuity 

The  continuity  equation  for  groundwater  flow  is 


V  •  q  =-S 


s  dt 


(3.6) 


in  which  S  ,  the  specific  storage,  is  defined  as 


where 


S  =  pg(a  +  n3) 


a 
6 
n 
P 


coefficient  of  soil  compressibility 
coefficient  of  compressibility  for  water 
porosity 
fluid  density 


(3.7) 


The  combination  of  equations  (3.2)  and  (3.6)  gives  the  governing 
equation  for  groundwater  flow 


_9_ 
9x 


' 
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(3.8) 


14 

For  steady  flow  in  homogeneous  and  isotropic  soil  this  reduces  to  the 
Laplace  equation 

i±,i±,ii-_^  (3.9) 

ax    8y    9z 

Well  Flow  in  Leaky  Aquifers 

In  dealing  with  flow  to  and  from  wells,  it  is  customary  to  denote 
the  head  gradient  in  terms  of  drawdown,  s,  which  is  defined  as 

s  =  (}>  -  (|)o  (3.10) 

where  4)q  is  the  head  above  the  leaky  layer  (see  Figure  3.1).  The  flow 
from  an  injection  well  causes  an  increase  of  head,  in  which  case  s  will 
be  called  "push  up"  instead  of  drawdown.  A  leaky  aquifer  is  one  which 
is  confined  by  a  layer  of  material  which  has  low  hydraulic  conductivity, 
but  not  so  low  that  it  can  be  considered  impervious.  The  governing 
equation  for  flow  in  a  leaky  aquifer  is 

The  leakage  coefficient,  B,  is  given  by 

B^  =  Kbb'/K'  (3.12) 


where 


K'  =  hydraulic  conductivity  of  the  leaky  layer 
b'  =  thickness  of  the  leaky  layer 
b  =  thickness  of  the  aquifer 
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The  Dispersion  Equation 

The  mixing  of  saltwater  with  freshwater  at  the  periphery  of  a 
freshwater  lens  is  described  by  the  convective  dispersion  equation 

9I+  -  q  •  Vc  =  V  .  D  .  Vc  (3.13) 

in  which  c  is  the  concentration  of  the  fluid  being  followed  and  D  is 
the  dispersion  tensor.  For  isotropic  porous  media  this  tensor  has  only 
two  independent  terms,  D^  and  D^,  corresponding  to  longitudinal  and 
transverse  dispersion  respectively,  relative  to  the  direction  of  the 
mean  flow.  The  values  of  these  coefficients  depend  not  only  on  the 
properties  of  the  porous  medium  but  also  on  the  flow  velocity  and  the 
properties  of  the  fluid  involved.  It  is  often  assumed  that  these 
coefficients  can  be  represented  by 

°£  =  ^m^^^q^/"  (3.14a) 

^t  =  °m  ^  «t  ^£/"  (3.14b) 

where 

D^  =  coefficient  of  molecular  diffusion  in  the  porous  medium 

a.   =  longitudinal  dispersivity 

a.  =  transverse  dispersivity 

q^  =  discharge  velocity  in  the  direction  of  mean  flow 

The  dispersivity  coefficients  a^  and  a^  are  properties  of  the 
porous  medium  only.  Their  magnitudes  are  related  to  the  intrinsic 
permeability  and  the  pore  size  distribution,  but  the  actual  values  must 
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be  found  experimentally.  The  longitudinal  dispersion  coefficient,  D  , 
is  usually  10  to  20  times  greater  than  D^  but  the  ratio  of  the  two 
depends  on  the  flow  velocity,  being  close  to  unity  when  the  flow  is 
very  slow  and  the  influence  of  molecular  diffusion  predominates. 

For  axi symmetric  radial  flow  of  a  tracer-like  substance  from  an 
injection  well,  equation  (3.13)  can  be  simplified  to 

at  *  ^ZTTbnr^  9?  -  «£  ^2^^  ^  (3.15) 

because  there  are  no  components  of  either  velocity  or  concentration 
gradient  except  in  the  radial  direction.  This  equation  has  been 
solved  numerically  for  the  case  of  steady  flow  from  the  well  and  con- 
tinuous tracer  injection  beginning  instantanously  (Hoopes  and  Harleman, 
1965;  and  Shamir  and  Harleman,  1966). 

Freshwater  injected  into  a  saline  aquifer  does  not  behave  as  a 
tracer  because  of  the  difference  in  density,  so  equation  (3.15)  does 
not  apply  even  though  the  flow  may  be  asixymmetric.  Because  vertical 
components  of  both  velocity  and  concentration  gradient  are  present, 
equation  (3.13)  must  be  used.  The  numerical  solution  of  this  equation 
will  be  the  subject  of  Chapter  7  of  this  report. 

Boundary  Conditions  at  an  Interface 
When  the  mixing  zone  separating  saltwater  and  freshwater  is 
relatively  thin,  it  is  often  convenient  to  neglect  dispersion  and 
assume  that  the  two  fluids  are  separated  by  a  sharp  interface.  If 
such  a  boundary  between  the  two  fluids  exists,  a  separate  governing 
equation  can  be  written  for  each  zone  and  the  two  flows  will  be  coupled 
by  their  common  boundary  condition  at  the  interface. 
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The  interface  is  a  material  surface  which  is  always  composed  of 
the  same  fluid  particles.  The  essential  features  of  such  a  surface 
are  that  both  pressure  and  the  velocity  component  normal  to  the 
interface  must  be  continuous  across  it.  For  a  stationary  interface 
the  normal  velocity  component  is  zero  so  the  boundary  condition  is 
concerned  with  pressure  only.  The  pressure  on  the  freshwater  side  of 
the  interface  is 

Pf  =  Yf  {<l>^  -   z.)  (3.16a) 

and  on  the  saltwater  side 

Pg  =  ^s  ^*s  "  ^-^  (3.16b) 

where  z^.  is  the  elevation  of  a  point  on  the  interface  and  the  sub- 
scripts f  and  s  stand  for  freshwater  and  saltwater,  respectively. 
Requiring  that  these  pressures  be  equal  results  in  an  expression  for 
the  boundary  condition  at  a  stationary  interface 


Y       Y 

'i  ^  a7  *s  -  a7  *f  (3.17) 

where 


AY  =  Y3  -  Y^. 

For  the  special  case  where  the  saltwater  is  at  rest  and  a  hydro- 
static pressure  distribution  prevails  in  the  freshwater  zone  this 
reduces  to  the  Ghyben-Herzberg  relationship 

^f 
Az  =  -^A4)^  (3.18a) 
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or 


^z  =  -  6A(f)^  (3.18b) 

where  6  =  y^/Ay  is  the  Ghyben-Herzberg  ratio.  This  states  that  a  change 
in  the  piezometric  head  of  the  freshwater  causes  a  change  in  the  eleva- 
tion of  the  interface  that  is  opposite  in  sign  and  greater  in  magnitude 
by  a  factor  of  6. 

If  the  interface  is  not  stationary  the  condition  of  continuity  of 
the  normal  velocity  component  must  also  be  used.  The  interfacial  surface 
can  then  be  described  by  some  function,  F,  such  that 

F  (x,  y,  z,  t)  =  0  (3.19) 

This  surface  is  convected  by  the  flow  so  that 

lr  +  ^VF=0  (3.20) 

where 

q  =  discharge  velocity  of  interfacial  particles 
n  =  porosity  of  the  aquifer 

In  terms  of  the  velocities  adjacent  to  the  interface  on  both  sides  of  it 
this  becomes 


n  g^+  q^  •  VF  =  0  (3.21a) 


for  the  freshwater,  and 


"  at  ^^s  •  ^^  =  0  (3.21b) 
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for  the  saltwater.     In  terms  of  the  elevation  of  a  point  on  the  inter- 
face,  z^,  the  function,  F,  can  be  written  as 

F{x,  y,  z,   t)   =  z  -  z.   (x,  y,   t)   =  0  (3.22) 

Substituting  this,   together  with  Darcy's  law,  into  equations  (3.21)   gives 

9Z, 


■"  TT  "  "^^^f  ■   ^^^  -  Zi)  =  0  (3.23a) 


8z. 


■n  —-  -  Kv^^   •   v(z  -  z.)   =  0  (3.23b) 


Performing  the  vector  operations  indicated  in  these  equations  results 
in  the  boundary  conditions  for  the  freshwater  zone  and  the  saltwater 
zone.  The  interface  boundary  condition  for  the  freshwater  zone,  then, 
is 


Yf  8<t'f     y     90      Yf     n  y. 


-K-^=0  (3.24a) 


and  for  the  saltwater  zone. 


Yf  9-1)^       y      9({>        Y<-       p       Yf 


K^=0  (3.24b) 
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Fluid  Properties 
Density 

Up  to  this  point  the  fluids  with  which  this  report  is  concerned 
have  been  referred  to  simply  as  freshwater  and  saltwater.  These  conven- 
ient  terms  are  used  only  in  a  relative  sense  and  the  actual  properties 
of  the  fluids  deserve  some  further  elaboration.  In  the  ensuing  analy- 
sis the  primary  difference  between  the  resident  and  injected  fluids 
will  be  their  density.  In  general,  this  is  a  property  that  depends 
on  pressure,  temperature  and  salt  concentration.  In  this  report, 
however,  the  change  in  density  due  to  compressibility  will  be  neglected 
because  it  is  so  small,  and  the  thermal  effects,  while  not  necessarily 
small,  will  be  eliminated  by  assuming  that  the  temperature  of  the 
injected  fluid  is  the  same  as  that  of  the  resident  saltwater. 

Figure  3.2  shows  the  relationship  of  density  to  temperature  for 
water  with  several  different  salt  concentrations.  With  the  exception 
of  the  3.5%  concentration  line,  this  figure  is  based  on  data  for 
solutions  of  pure  sodium  chloride  in  water.  The  3.5%  concentration 
line  shows  the  temperature  versus  density  relation  for  seawater  having 
a  dissolved  solids  concentration  of  3.5%.  Saline  groundwater  is 
commonly  associated  with  the  intrusion  of  seawater  into  coastal  aquifers. 
However,  it  is  also  found  in  many  locations  which  have  no  apparent 
geological  association  with  the  oceans  and  where  its  concentration  may 
be  much  greater  or  less  than  that  of  seawater. 

Figure  3.3  shows  the  relationship  of  density  to  sodium  chloride 
concentration.  It  is  worthwhile  to  note  that  in  this  range  the 
relationship  is  apparently  linear. 
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FIGURE  3.2.  TEMPERATURE  VS.  DENSITY  FOR  WATER 
WITH  VARIOUS  SALT  CONCENTRATIONS 
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Viscosity 

It  was  shown  in  equation  (3.4)  that  the  hydraulic  conductivity 
of  an  aquifer  is  partially  determined  by  the  kinematic  viscosity  of  its 
water.  Viscosity  is  a  function  of  both  temperature  and  salt  concentra- 
tion but  within  the  range  of  present  interest  temperature  is  the  more 
important  factor  of  the  two.  The  temperature  dependence  of  kinematic 
viscosity  for  both  pure  water  and  seawater  is  shown  in  Figure  3.4.  It 
is  interesting  to  note  that  the  kinematic  viscosity  of  pure  sodium 
chloride  solutions,  as  given  by  Kaufmann  (1960)  for  concentrations  up 
to  4%,  differ  so  little  from  the  kinematic  viscosity  of  pure  water  that 
they  could  not  be  distinguished  from  it  if  plotted  at  the  scale  of 
Figure  3.4.  Since  differences  in  temperature  are  not  being  considered 
in  this  study  and  the  influence  of  salt  concentration  on  viscosity  is 
relatively  weak,  viscosity  differences  will  be  neglected  also. 
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FRESHWATER  AND  SALTWATER 


CHAPTER  4 
THE  HELE-SHAW  MODEL 

General 


Valuable  insights  into  the  flow  phenomena  associated  with  the 
injection  of  fresh  water  into  saline  groundwater  can  be  obtained  by 
studying  certain  simplified  aspects  of  the  problem.  This  chapter  will 
describe  the  use  of  a  physical  model  to  analyze  the  movement  of  the 
interface  between  the  two  fluids  in  the  absence  of  dispersion.  The 
Hele-Shaw  analog  was  chosen  for  this  work  because  it  is  simple  and 
inexpensive  to  construct  and  the  test  results  are  easily  observable 
without  elaborate  instrumentation.  The  Hele-Shaw  analog  consists  of 
one  or  more  viscous  liquids  flowing  in  the  gap  between  two  closely 
spaced  flat  plates.  Its  usefulness  in  groundwater  studies  is  based 
on  the  fact  that  the  Navier-Stokes  equations  governing  the  viscous 
flow  in  this  gap  can  be  reduced  to  a  form  similar  to  Darcy's  law  in 
two  dimensions. 

The  Viscous  Flow  Analogy 
If  u,  V,  and  w  are  the  velocity  components  in  the  x,  y,  and  z 
directions,  respectively  then  the  Navier-Stokes  equations  for  incom- 
pressible flow  in  a  cartesian  coordinate  system  are 
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i  ^"f  ^vf  -'H  =^z-^l^^vv2w       (4.10 


where 


b^,  b  and  b^  =  components  of  acceleration  due  to  body  forces  in 
the  X,  y,  and  z  directions,  respectively 

p  =  density  of  the  fluid 

V  =  kinematic  viscosity  of  the  fluid 

P  =  pressure 
Applying  these  equations  to  the  flow  between  the  plates  shown  in 
Figure  4.1,  it  can  be  assumed  that  because  of  the  low  Reynolds  number 
(Re  <  1  for  this  model)  the  flow  is  laminar  and  there  is  no  velocity 
component,  v,  in  the  direction  normal  to  the  plates.  Note  also  that, 
since  the  velocity  at  the  surface  of  the  plates  is  zero,  the  velocity 
gradients  normal  to  the  plates  are  \/ery   large  compared  to  those  in  the 
X  and  z  directions  so  that  the  latter  can  be  neglected.  Finally,  the 
only  body  force  acting  on  the  fluid  is  gravity,  so  b  =  -  g  and 
b  =  b  =  0.  Incorporating  these  observations  into  the  Navier-Stokes 
equations  and  multiplying  by  p  results  in 

9P  _   9^u  //,  o  \ 

33^  -  y-2  (4.2a) 


f  =  0  (4.2b) 


Y  +  -^  =  y  — ?  (4.2c) 

3y 


where  y  =  P9  is  the  unit  weight  of  the  fluid  and  y  =  pv  is  its  dynamic 
viscosity.  Recalling  that  the  piezometric  head,  4),  is  defined  as 
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*  =  z  +  P/y  (4.3) 

it  follows  that  equations  (4.2)  when  expressed  in  terms  of  (j)  and  divided 
by  Y  become 

^    -    77T  (4.4a) 

dy 

t  =  0  (4.4b) 

From  the  second  of  these  equations  it  is  seen  that  at  any  point  (x,  z) 
in  the  plane  of  the  model,  the  piezometric  head  is  not  a  function  of  y. 
Therefore,  the  first  and  third  of  these  equations  can  be  integrated 
with  respect  to  y  giving 

^  9x    Y  3y  (4.5a) 


^  9z  ~  7  ay  ^  ^  (4-5b) 

The  constant  of  integration,  C,  is  found  to  be  zero  by  observing  that 

at  the  center  of  the  flow  region,  where  y  =  0  both  -^  and  ~  are  also 

dy    dy 

zero,  so  that 


•y  3x    y  dy  ^4.6aj 


^  9z    Y  3y  (4-^^^ 


Integrating  these  equations  again  with  respect  to  y  gives 


-1 
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2 
V  f    Y  ^  -^  ^1  (4.7b) 

Here,  the  constant  of  integration,  Cp  is  found  by  requiring  that,  at 
the  surface  of  the  plates,  where  y  =  ±  d/2,  u  =  w  =  0.  With  these 
boundary  conditions  the  expressions  for  the  velocity  distribution 
across  the  gap  are  found  to  be 

u  =  i  (y^  -  ^)  t  (4.8a) 

where  d  is  the  thickness  of  the  gap  between  the  plates.     Now,   in  order 
to  find  a  velocity  analogous  to  Darcy's  discharge  velocity  for  flow 
through  porous  media,  equations   (4.8)  must  be  averaged  across  the 
thickness  of  the  gap  as  follows 


/.d/2 


d/2  2 


-/ 


d/2 
d/2 


-d/2 

If  the  constant  coefficients   in  equations    (4.9)   are  used  to  define  an 
effective  hydraulic  conductivity,  K,^,  for  the  model 

.2 
'<m  =  ^J  (4.10) 
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then  the  governing  equations  for  flow  between  the  plates  bee 


ome 


%--^^^  (4.11a) 


^z  =  -  ^m|f  (4.11b) 


in  which  the  resemblance  to  the  two  dimensional  form  of  Darcy's  law  is 
obvious.  In  further  analogy  to  groundwater  flow,  an  intrinsic  permea- 
bility for  the  Hele-Shaw  model,  k^,  can  be  defined  as 


so  that 


km=d2/l2  (4J2) 


•^m  =  '<m  J  (4.13) 


The  continuity  equation  for  the  flow  between  the  plates  is 


9^  ^  91"  0  (4.14) 


since  there  is  no  velocity  in  the  y  direction.  Averaging  these  veloc- 
ities across  the  gap,  as  before,  results  in 

ir  (^m  If)  *  i  (><- If)  =  0  (4.15) 

Which  has  the  same  form  as  the  governing  equation  for  two  dimensional 
flow  in  an  incompressible  porous  medium.  The  model  hydraulic  conductiv- 
ity, K^,  is  a  constant  depending  only  on  the  spacing,  d,  of  the  plates 
and  the  unit  weight  and  viscosity  of  the  fluid.  Consequently,  for 
parallel  plates  and  homogeneous  fluids  the  flow  is  homogeneous  and 
isotropic,  so  that  equation  (4.15)  reduces  to  the  Laplace  equation  in 
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two  dimensions 


-^  +  -^  =  0  (4.16) 


Since  K^  is  a  function  of  physical  parameters  which  are  relatively 
easy  to  measure  and  control,  the  Hele-Shaw  model  is  a  very  convenient 
tool  for  the  scale  modelling  of  two  dimensional  groundwater  flow. 

Model  Scaling 

In  the  preceding  section  it  was  shown  that  two  dimensional  ground- 
water flow  can  be  modeled  by  the  Hele-Shaw  analog  because  the  governing 
equations  for  the  two  cases  are  of  the  same  form.  If  quantitative 
results  are  to  be  obtained  from  the  model,  however,  the  scaling  ratios 
between  model  and  prototype  must  be  chosen  so  that  the  governing 
equations  are  actually  identical.  In  order  to  achieve  this,  the  model 
laws  which  fit  the  Hele-Shaw  model  to  the  prototype  case  of  groundwater 
flow  must  be  developed. 

Since  the  present  case  involves  the  modeling  of  two  immiscible 
liquids  separated  by  an  interface  the  equations  on  which  the  model 
laws  depend  are  Darcy's  law,  the  continuity  equation  and  the  boundary 
condition  for  interface  flow.  Letting  the  subscript  m  denote  model 
parameters  and  p  denote  those  in  the  prototype,  Darcy's  law  for  the 
model  is 


m 
m 
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(  f  for  the  lighter  fluid 
where  the  superscript        i  =      < 

(s  for  the  heavier  fluid 


The  corresponding  equations  for  the  prototype  are 


d<t>] 


q'     =-Kp  9/  (4.18a) 

xp  ^  ''^p 

i 

'zp  "p   3z^ 


9(f) 

q'n    =    -K   ^  (4.18b) 


Let  the  scaling  ratios  between  model   parameters  and  prototype  parameters 
be  denoted  by  the  subscript  r,  so  that 


K,=  VKp  {4.19a) 


q^  =  V^p  (4.19b) 


\  =  VXp  (4.19c) 


etc. 


Then,  substituting  these  relationships  into  equations  (4.18)  gives 


Kn,  x^  a*' 


<Xr----i-  (^-^Oa) 


^  *r  '\ 


K     1      90^ 
.1  /„i      m  r  ^m 


W^r=---^—  (4.20b) 

r  ^r   r 

If  these  are  to  be  identical  to  equations  (4.17)  then  it  must  be 

ture  that  ,.  i 

i    "^r^r 
\r  "  ^  (4.21a) 
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Equations  (4.21)  provide  one  of  the  model  laws  which  allow  results 
from  the  model  to  be  transferred  to  the  prototype. 

The  continuity  equations  for  isotropic  flow  in  the  model  and 
prototype  are 

-^    ^    -^    -0  (4.22) 

^^m     K 


and 


2  2 


(4.23) 


Substituting  the  appropriate  parameter  ratios  into  equation  (4.23)  and 
requiring  that  it  be  identical  to  equation  (4.22)  results  in 

\  =  ^r  (4.24) 

The  third  similitude  requirement  is  that  the  boundary  conditions 
at  the  moving  interface  be  given  by  identical  equations  in  both  model 
and  prototype.  These  interface  equations  were  given  in  Chapter  3.  For 
the  model  they  are 

-  K^3^  =0  (4.25a) 

m 
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and 


Y^  9/     y'  acf)'     y'     2     Y^    . 


i 


m   9z 
m 


=  0  (4.25b) 


The  corresponding  equations  for  the  prototype  are  the  same  as  equations 
(4.25)  except  that  the  m  subscripts  are  replaced  by  p.  When  the  param- 
eter ratios  are  substituted  into  the  prototype  equations  and  equality 
is  required  with  the  model  equations,  16  groupingsof  parameter  ratios 
are  generated,  which  must  all  be  equal  to  unity.  Many  of  these  group- 
ings are  repetitive  and,  since  they  provide  no  new  information,  they 
will  not  be  listed  here.  The  groupings  which  are  useful  are 


1  AY^  t     1  AY^  t     1  AY„  xl  1  z     1  2 

n,  Y,  *,    n   y^c^^     K   yj  ^\  Y.^  ^\  K^  ^^ 


The  equality  of  the  first  and  third  terms  of  this  equation  gives  the 
time  scale  for  the  model 

"r  \ 
V  =  ^^  (4.27) 

r  ^r 

Equality  of  the  first,  second,  and  fourth  terms  of  equation  (4.26) 
together  with  equation  (4.27)  results  in 


f  f    s  s 

Y^  (t>p  ^  '^y^^     =   Ay^  Zp  (4.28) 
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Finally,  the  equality  between  the  last  two  groupings  in  equation  (4.26) 
requires  that 

which,  in  turn,  means  that 

f     s 

Yp  =  Yy,  (4.30) 

The  complete  list  of  similitude  requirements  for  matching  the 
prototype  to  the  model  is  given  in  Table  (4.1). 

Model  Construction 
The  design  of  the  Hele-Shaw  model  involved  the  following  con- 
sideration: 

1.  The  plates  must  be  transparent,  relatively  inflexible  and 
demonstrably  flat. 

2.  Provision  must  be  made  for  the  controlled  introduction  of  the 
viscous  fluids  and  the  release  of  entrapped  air. 

3.  The  model  must  be  easily  disassembled  for  cleaning  between 
tests. 

4.  It  must  be  possible  to  rotate  the  model  from  horizontal  to 
vertical  so  that  the  initial  shape  of  the  interface  can  be 
controlled. 

The  availability  of  materials  and  the  requirement  for  transparency 
dictated  that  the  material  to  be  used  for  the  flat  plates  be  either 
1/4  inch  plate  glass  or  1/2  plexiglass.  Since  1/2  plexiglass  is  more 
than  2.5  times  more  flexible  than  1/4  inch  plate  glass,  it  was  decided 
to  use  the  latter. 
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TABLE  4.1 

SUMMARY  OF  SIMILITUDE  REQUIREMENTS  FOR  VERTICAL 
HELE-SHAW  MODEL  WITH  TWO  IMMISCIBLE  FLUIDS 


Velocity  Ratios; 


i    ^  *r  i     K  < 

q    =  ^r-^     and     q'   ^  -L^ 

xr     \  zr     ^r 


(f,   for  the  lighter  fluid 
where     i  =  <[ 

(s,   for  the  heavier  fluid 
Geometric  Ratios: 


\  =  ^r 


Time  Ratio: 

2 

•"  '  K    / 


r  ^r 


Head  and  Density  Ratios: 


f   A.^   -      s  ,  s 

\  *r  -  ^r  *r  =  '^^r  ^r 


and 


f     s 

'r     'r 
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It  was  originally  thought  that  all  plate  glass  was  \iery   flat  and 
that  its  flatness  must  be  controlled  to  very   close  tolerances  during 
its  manufacture.  Conversations  with  production  engineers  at  several 
glass  factories  revealed,  however,  that  most  emphasis  is  placed  on  the 
uniformity  of  thickness  of  the  glass  and  not  on  its  flatness.  Essen- 
tially, all  standard  plate  glass  is  now  produced  by  the  flotation 
process.  This  gives  it  a  very  smooth  surface  and  uniform  thickness  and 
it  is  simply  assumed  to  be  flat  enough  for  most  purposes.  Numerous 
pieces  of  commercially  available  plate  glass  were  optically  tested  for 
flatness,  and  it  was  found  that  some  pieces  are  much  flatter  than 
others.  The  tests  were  run  by  placing  the  carefully  cleaned  surfaces 
of  two  glass  plates  together  and  illuminating  them  with  monochromatic 
light  from  a  sodium  vapor  lamp.  The  uniformity  of  the  interference 
fringes  produced  gives  a  qualitative  indication  of  the  flatness  of  the 
two  glass  panes.  Three  pieces  of  special  plate  glass  were  obtained 
which  were  manufactured  by  the  traditional  process  of  mechanical  grind- 
ing and  polishing.  The  manufacturer  claimed  that  this  glass  was  flat 
to  within  two  interference  fringes  per  inch,  where  one  interference 
fringe  represents  a  deviation  from  the  mean  surface  of  approximately 
2.75  X  10   millimeters.  When  these  pieces  of  plate  glass  were  observed 
under  the  sodium  vapor  lamp,  their  superior  flatness  was  evident. 
Selected  pieces  of  the  commercially  available  glass,  however,  were 
nearly  as  good.  These  were  used  in  the  construction  of  the  Hele-Shaw 
model . 

In  order  to  provide  for  the  controlled  infusion  of  the  viscous 
liquids  into  the  gap  between  the  plates,  it  was  necessary  to  glue  a 
section  of  machined  plexiglass  to  each  end  of  both  glass  plates.  In 
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this  way,  an  infusion  reservoir  was  created  on  either  end  of  the  model 
so  that  the  liquids  could  be  injected  uniformly  into  the  test  section. 
This  is  shown  in  Figure  4.2. 

The  gap  between  the  two  glass  plates  was  maintained  at  a  uniform 
thickness  by  a  series  of  brass  spacers  1/16  of  an  inch  (1.59  mm)  thick 
placed  at  intervals  around  the  edges  of  the  model.  The  plates  were 
held  in  position  by  C-clamps  which  pressed  them  tightly  against  the 
spacers.  The  model  was  sealed  by  a  length  of  surgical  rubber  tubing 
which  was  placed  around  the  edges,  just  inside  of  the  spacers,  and 
crushed  flat  by  the  action  of  the  C-clamps.  The  wall  thickness  of  this 
tubing  was  approximately  1/32  of  an  inch  so  that  it  could  be  crushed 
between  the  plates  without  exerting  much  force  against  the  glass. 

The  assembled  model,  mounted  on  its  stand  is  shown  in  Figure  4.3. 
The  length  of  the  test  section  is  76  centimeters  and  its  height  is  38 
centimeters.  With  the  rubber  seal  in  place,  however,  the  actual  height 
available  for  flow  is  35.6  centimeters. 

The  viscous  liquids  are  supplied  to  the  infusion  reservoirs  on 
either  end  of  the  model  through  plastic  tubes.  The  plexiglass  end 
pieces  on  the  back  plate  are  each  supplied  with  four  fittings  for  1/2 
inch  (1.27  cm)  plastic  tubing.  The  heavier  fluid  is  conducted  to  the 
model  through  this  tubing  from  constant  head  tanks  mounted  on  a  ring 
stand.  This  is  shown  in  Figure  4.4.  On  the  front  plate,  the  end 
pieces  are  provided  with  air  vents  on  one  end,  and  four  nipples  for 
the  connection  of  1/4  inch  (0.64  cm)  plastic  tubing  on  the  other. 
The  lighter  viscous  fluid  is  injected  through  this  tubing  into  one  of 
the  infusion  reservoirs  by  a  positive  displacement  syringe  pump.  This 
setup  is  illustrated  in  Figure  4.5. 


40 


LU 
Q 

a: 


LU 
LU 

a: 


I — I 
CO 


_i 

Q. 


LU 

=3 
Ci3 


41 


■/:Ah 


m 


< 

I— 


o 
O 

ZSl 

Q 

UJ 

-J 
in 

CO 
CO 

CO 


a: 

ID 


42 


Q 

■a: 

UJ 


■a: 

I— 

00 


o 


Q 
O 


LlJ 
C3 


O 
GQ 


LU 


43 


44 


Properties  of  the  Viscous  Fluids 
Dow  Corning  Corporation  Series  200  silicone  oil  was  used  as  the 
heavier  fluid  in  the  model.  This  is  a  clear  dimethyl  siloxane  which 
is  characterized  by  oxidation  resistance  and  low  vapor  pressure.  Its 
specific  gravity  is  0.977  and  its  kinematic  viscosity  is  nominally  500 
centistokes  at  25°C.  For  the  lighter  fluid  in  the  model,  a  mixture  of 
90  SAE  and  30  SAE  non-detergent  motor  oils  was  used.  This  mixture  was 
adjusted  until  its  kinematic  viscosity  was  the  same  as  that  of  the 
silicone  oil  at  a  temperature  of  24.5°C.  The  temperature  versus 
viscosity  relationships  for  the  two  fluids  in  the  neighborhood  of  25°C 
is  shown  in  Figure  4.6.  This  data  was  obtained  with  a  Cannon-Zietfuchs 
viscometer  with  calibration  traceable  to  the  National  Bureau  of  Stan- 
dards. The  specific  gravity  of  the  light  oil  mixture  was  0.890  at 
25°C. 

These  two  fluids  are  not  readily  miscible,  yet  the  surface  tension 
effects  between  the  two  are  not  so  great  that  capillary  action  at  the 
interface  is  noticeable. 

Operation  of  the  Model 
The  Hele-Shaw  model  is  a  versatile  analytical  tool  and  its  mode 
of  operation  depends  somewhat  on  the  type  of  test  being  run.  In  the 
present  case  the  object  of  the  tests  was  to  observe  the  shape  and  rate 
of  rotation  of  the  initially  vertical  interface  between  the  two  fluids 
of  different  density.  The  idea  was  that  the  behavior  of  such  an 
interface  would  be  similar  to  that  of  the  interface  formed  by  the 
injection  of  freshwater  from  a  well  into  a  confined  saline  aquifer, 
neglecting  dispersion.  Of  course,  the  two  processes  are  not  the  same 
because  the  flow  away  from  a  well  is  radial  while  the  Hele-Shaw  model 
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is  limited  to  two-dimensional  flow  in  cartesian  coordinates.  Direct 
application  of  the  Hele-Shaw  results  to  a  prototype  well  situation  is, 
therefore,  not  possible.  With  increasing  distance  from  the  well, 
however,  the  two  flows  become  more  and  more  similar  so  that  the  model 
tests  are  useful  for  the  qualitative  description  of  interface  movement. 
The  results  of  the  model  tests  are  also  useful  for  comparison  with 
numerical  models  which  can  be  made  to  work  in  either  radial  or 
cartesian  coordinate  systems. 

In  setting  up  the  model  for  a  test,  it  is  necessary  first  to  be 
sure  that  all  parts  of  it  are  clean  and  dry.  This  can  only  be  done  by 
completely  disassembling  it  and  washing  each  part.  A  commercial  pre- 
paration of  trichloroethane  cleaning  solution  was  found  useful  in 
washing  the  oils  from  the  glass  and  plastic  parts.  After  the  model  had 
been  cleaned  and  reassembled,  with  the  plates  rotated  to  the  horizontal 
position,  the  heavier  oil  was  introduced  at  the  downstream  end;  the 
upstream  end  being  the  one  into  which  the  lighter  oil  would  be  injected. 
During  this  process  the  upstream  end  was  raised  a  few  centimeters  so 
that  all  of  the  air  in  the  test  section  would  be  expelled.  Air 
entrapped  in  the  downstream  infusion  reservoir  was  released  by  loosen- 
ing the  brass  plugs  provided  for  this  purpose.  The  introduction  of 
silicone  oil  was  continued  until  the  oil  reached  the  end  of  the  glass 
and  was  just  ready  to  run  into  the  upstream  infusion  reservoir.  At 
this  point,  the  flow  was  stopped  and  the  oil  front  was  held  at  the 
edge  of  the  glass  by  proper  adjustment  of  the  constant  head  supply 
tank.  The  1/2  inch  (1.27  cm)  plastic  tubes  on  the  bottom  side  of  the 
upstream  end  were  then  clamped  off  tightly  with  hose  clamps.  Three 
of  the  1/4  inch  (0.635  cm)  tubes  from  the  syringe  pump  were  then 
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hooked  up  and  the  upstream  infusion  reservoir  was  filled  with  motor 
oil.  The  fourth  nipple  was  left  open  for  the  escape  of  air.  After 
all  of  the  air  had  escaped  and  the  model  was  entirely  filled  with  the 
two  viscous  fluids,  the  fourth  tube  was  attached  and  the  interface 
between  the  two  oils  was  driven  into  the  model  to  the  starting  point 
for  the  test.  It  was  found  helpful  to  let  the  interface  sit  in  this 
position  for  several  hours  so  that  any  irregularities  in  it  would  be 
straightened  out  by  gravity.  The  elevated  end  of  the  model  was  then 
lowered  to  a  level  position  and  the  test  begun. 

The  test  was  started  by  simultaneously  rotating  the  model  to  the 
vertical  position,  turning  on  the  syringe  pump  and  starting  the  timer. 
The  movement  of  the  interface  could  then  be  photographed  until  the 
downstream  end  of  the  model  was  reached. 

Several  variations  of  this  procedure  were  used.  The  first  test 
was  run  with  the  model  left  in  its  horizontal  position  so  that  the 
straightness  of  the  advancing  interface  could  be  checked.  It  was  found 
that,  with  higher  rates  of  injection  of  the  lighter  fluid,  the  interface 
moved  faster  in  the  center  of  the  model.  This  indicated  that  the 
higher  pressure  required  to  drive  the  viscous  fluids  across  the  model 
could  cause  the  glass  plates  to  bend,  changing  the  thickness  of  the 
gap.  Figures  4.7  and  4.8  show  the  shape  of  the  interface  when 
advancing  at  rates  of  3.65  cm/mi n  and  7.28  cm/mi n,  respectively.  Since 
the  straightness  of  the  interface  at  the  lower  speed  was  considered 
acceptable  most  of  the  remaining  tests  were  run  at  this  rate  of 
injection. 

Figures  4.9  and  4.10  show  a  test  which  was  run  at  a  higher  injec- 
tion rate  with  the  model  in  the  vertical  position.  The  bowing  out  of 
the  glass  plates  caused  the  center  section  of  the  interface  to  advance 
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more  rapidly  than  the  ends.  Rotation  of  the  upper  part  of  the  interface 
past  the  vertical  results  in  instability,  which  is  clearly  evident  here. 

The  sequence  of  Figures  4.11  through  4.16  shows  a  test  that  was 
run  with  an  injection  rate  of  20.44  cubic  centimeters  per  minute. 
This  is  the  same  injection  rate  that  resulted  in  the  relatively  straight 
interface  shown  in  Figure  4.7  when  the  model  was  horizontal. 

The  sequence  of  Figures  4.17  through  4.23  shows  a  test  in  which 
only  the  lower  half  of  the  model  was  used.  This  was  accomplished  by 
inserting  a  barrier  composed  of  rubber  splicing  compound  and  silicon 
rubber  sealant  between  the  glass  plates  to  reduce  the  width  of  the  test 
section  by  half.  The  injection  rate  used  here  was  8.16  cubic  centi- 
meters per  minute. 

Analysis  of  Test  Results 
The  rotation  of  an  initially  vertical  interface  between  two  fluids 
of  different  density  has  been  studied  in  several  previous  investigations. 
Gardner,  Downie  and  Kendall  (1962)  proposed  an  equation  to  predict  the 
rate  of  interface  rotation  which,  for  an  isotropic  aquifer  and  fluids  of 
equal  viscosity,  has  the  form 


/X,^  _  16  (t/t')^  ,,   „., 

^F    ~Xin/t'  ^^-31) 


where 

X  =  projection  of  the  interface  on  the  horizontal 

b  =  thickness  of  the  aquifer 

t  =  time 

f    =  4/3  ^ 

p  =  average  density  of  the  two  fluids. 
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Equation  (4.31)  was  constructed  intuitively  by  joining  the  solutions  to 
the  extreme  cases  of  nearly  vertical  and  nearly  horizontal  interfaces. 
The  solutions  to  these  extreme  cases  were  based  on  the  assumptions  of 
no  vertical  flow  and  no  curvature  of  the  interface. 

Another  equation  for  the  movement  of  an  initially  vertical  inter- 
face was  given  by  Harleman  and  Rumer  (1962).  Once  again,  assuming  a 
straight  interface  and  no  vertical  flow,  they  obtained 


X  2      t 
(^)  =  8/3  |.  (4.32) 


This  equation  is  of  the  same  form  as  the  one  for  nearly  horizontal 
interfaces  from  which  equation  (4.31)  was  derived. 

Figure  (4.24)  shows  the  comparison  of  equations  (4.31)  and  (4.32) 
with  data  taken  from  three  runs  of  the  vertical  Hele-Shaw  model.  Runs 
number  two  and  three  are  shown  in  Figures  4.11  through  4.16  and  4.17 
through  4.23,  respectively.  Run  number  one  was  performed  in  the  same 
manner  as  run  number  two  except  that  the  injection  rate  was  twice  as 
high.  The  third  run  gave  the  best  fit  to  equation  (4.31).  Runs  number 
one  and  two,  while  fairly  consistent  between  themselves,  do  not  fit 
either  equation  very  well.  Since  the  main  difference  between  these  two 
runs  and  the  third  run  is  in  the  depth  of  the  flow  region,  it  is  most 
logical  to  suspect  that  it  is  caused  by  the  existence  of  vertical 
velocity  components  which  were  ignored  in  the  derivation  of  the 
equations.  The  greater  curvature  of  the  interface  which  is  apparent  in 
the  first  two  runs  tends  to  support  this  view. 
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FIGURE  4.24.  COMPARISON  OF  EXPERIMENTAL  RESULTS  WITH 
EQUATIONS  FOR  INTERFACE  ROTATION 


CHAPTER  5 
SIMPLE  CASES  INVOLVING  STEADY  FLOW 

General 


It  was  shown  in  Chapter  3  that  the  general  problem  of  locating  a 
moving  interface  between  two  fluids  of  different  density  is  very  complex 
even  when  the  standard  simplifying  assumptions  involving  homogeneity, 
isotropy  and  immiscibility  are  used.  It  requires  the  solution  of  two 
coupled  flow  problems,  one  in  the  freshwater  and  one  in  the  saltwater, 
joined  by  a  nonlinear  boundary  condition  on  the  interface,  of  which  the 
position  is  unknown.  The  only  methods  presently  available  for  getting 
a  solution  to  such  a  problem  are  by  physical  and  numerical  modeling. 
However,  if  the  problem  is  limited  to  the  special  case  of  a  stationary 
freshwater  lens  in  a  static  saline  aquifer,  then  it  becomes  much  more 
tractable.  In  this  special  case,  it  is  only  necessary  to  solve  the 
continuity  and  motion  equations  in  the  freshwater  region  and  the  location 
of  the  interface,  while  still  not  initially  known,  can  be  treated  by  the 
Ghyben-Herzberg  principle.  This  makes  it  possible  to  formulate  some 
problems  which  can  be  solved  analytically. 

Steady  flow  in  a  freshwater  lens  is  not  possible  unless  the  aquifer 
possesses  a  constant  head  boundary  which  permits  freshwater  to  flow 
out  of  the  lens  at  the  same  rate  that  it  is  being  injected.  One  such 
constant  head  boundary  is  found  at  the  sea  coast  where  injection  wells 
may  be  used  to  create  freshwater  lenses  in  aquifers  affected  by  saltwater 
intrusion.  Another  case  is  that  of  injection  into  a  leaky  saline  aquifer 
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where  the  leaky  layer  serves  as  a  boundary  which  permits  efflux  of 
freshwater  from  the  lens.  Even  in  these  cases,  steady  flow  may  only 
occur  during  a  small  part  of  the  history  of  the  lens.  It  represents  a 
condition  of  equilibrium  in  which  buoyant  forces  are  counteracted  by  the 
loss  of  energy  due  to  the  flow  resistance  of  the  porous  medium.  It  is 
useful  to  study  this  condition  because  it  indicates  the  limiting  size 
of  lens  that  can  be  achieved  under  the  given  flow  conditions. 

It  should  also  be  noted  that,  in  the  operation  of  a  storage  lens, 
certain  advantages  can  be  gained  by  maintaining  the  lens  in  a  steady 
state  through  continuous  freshwater  injection.  The  first  is  that 
continuous  injection  is  necessary  to  preserve  the  potential  energy  of 
the  lens  so  that  the  maximum  amount  of  freshwater  can  be  withdrawn  when 
it  is  needed.  Note  that  as  soon  as  injection  is  stopped,  the  lens 
begins  to  decay  into  a  thin  layer  of  freshwater  at  the  top  of  the  aqui- 
fer which  could  not  be  selectively  retrieved  even  if  fresh  and  saltwater 
did  not  mix.  The  second  advantage  is  that  the  continuous  injection  of 
freshwater  to  the  lens  tends  to  wash  out  the  saltwater  which  had  re- 
mained in  the  pores  and  adhered  to  the  particles  of  the  porous  medium 
after  the  passage  of  the  mixing  zone  as  the  lens  was  being  formed.  In 
this  way,  the  salinity  of  the  freshwater  and  the  thickness  of  the 
mixing  zone  are  reduced. 

Freshwater  Injection  in  a  Coastal  Aquifer 
Saltwater  Intrusion 

It  is  commonly  observed  in  coastal  aquifers  which  discharge  to  the 
sea  that  a  wedge  of  nearly  static  saltwater  extends  inland  from  the 
coast  underneath  the  flowing  freshwater.  Actually,  it  was  the  study  of 
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this  phenomenon  which  led  W.  Badon  Ghyben  and  A.  Herzberg  to  develop 
their  hydrostatic  theory  for  interfaces.  Numerous  subsequent  investi- 
gators have  brought  forth  more  refined  methods  for  analyzing  saltwater 
wedges  (Glover,  1959;  Henry,  1959;  Bear  and  Dagan,  1964;  Van  Der  Veer, 
1977)  and  have  found  the  Ghyben-Herzberg  results  to  be  quite  accurate 
except  for  that  part  of  the  interface  closest  to  the  sea. 

The  application  of  this  principle,  combined  with  the  Dupuit  assump- 
tion, to  a  confined  aquifer  of  constant  thickness  gives  the  following 
relationship  between  h  and  x. 

X  =  Kh2/(26Q^)  (5.1) 

in  which  the  symbols  are  as  shown  in  Figure  5.1.  The  position  of  the 
toe  of  the  interface  is  given  by 

X   =^  (<,9\ 

h       26Q  (5.2) 

^x 

Injection  Into  a  Coastal  Aquifer 

If  an  injection  well  pumps  freshwater  at  a  constant  rate  into  the 
coastal  aquifer  pictured  in  Figure  5.2  then  the  equilibrium  position  of 
the  saline  wedge  will  be  modified  and  a  steady  freshwater  lens  will  be 
formed  in  what  was  formerly  a  saline  part  of  the  aquifer.  This  is  a 
case  of  injection  into  an  aquifer  in  which  only  the  freshwater  is  flow- 
ing, so  that  the  Ghyben-Herzberg  principles  can  still  be  used.  The  flow 
pattern  around  a  well  in  a  uniformly  flowing  aquifer  is  customarily 
studied  by  the  methods  of  potential  theory  in  which  the  velocity  poten- 
tials for  source  flow  and  uniform  flow  are  superimposed.  However,  in 
the  presence  of  the  saltwater  wedge,  even  though  the  vertical  velocity 
components  are  neglected,  flow  is  not  two  dimensional  because  of  the 
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changing  thickness  of  the  flow  region.  Strack  (1976)  has  found  a 
modified  form  of  the  vertically  integrated  potentials  originally  pro- 
posed by  Girinskii  (see  Bear,  1972;  p.  157)  which  can  account  for  the 
presence  of  the  interface.  If  the  thickness  of  the  flow  region  is  a 
linear  function  of  the  piezometric  head,  {}>,  i.e., 

h  =  atf)  +  e  (5.3) 

then  Strack  gives  the  potential  for  the  interface  zone  as 

<D.  =  l/2aK(4.  +  6/a)^  (5.4) 

For  the  inland  part  of  the  aquifer  where  there  is  no  interface,  the 
potential  is  defined  by 

f  =  Kb(j)  +  c  (5.5) 

where  the  constant,  c,  is  used  to  make  the  potential  function  continuous 
at  the  toe  of  the  interface.  It  is  easily  verified  from  equations  (5.4) 
and  (5.5)  that  the  discharge  components  in  the  x  and  y  directions  are 
obtained  from  $  in  the  usual  way. 

Q^  =  -9$/8x    and    Q^  =  -d'P/dy  (5.6) 

in  which  Q^  and  Q^  represent  the  discharge  through  a  unit  width  and  the 
total  thickness  of  the  aquifer  in  the  x  and  y  directions,  respectively. 
From  equations  (5.6)  it  follows  that  the  continuity  condition  for  the 
confined  coastal  aquifer  is  given  by  the  Laplace  equation 

9^<I>/8x^  +  8^$/8y^  =  0  (5.7) 

and  $  is  a  harmonic  function  of  x  and  y.     The  effect  of  this   is  that 
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through  use  of  Girinskii's  potential,  variations  in  the  vertical  direc- 
tion have  been  integrated  out  of  the  problem.     The  resulting  potential 
function  can  be  manipulated,   using  the  theory  of  complex  variables  to 
describe  a  two-dimensional   flow  field  with  the  desired  boundary  condi- 
tions.    At  the  same  time,  the  numerical    value  of  the  potential   function 
at  any  point  gives  the  depth  to  the  interface  at  that  point  through 
the  Ghyben-Herzberg  relation  and  thereby  supplies  three-dimensional 
information  about  the  freshwater  lens.     Before  this  can  be  done,  however, 
it  is  necessary  to  determine  the  constants  a  and  6  in  equation  (5.4)  and 
c  in  equation  (5.5) . 

Confined  coastal   aquifer.      From  the  geometry  of  the  saltwater  wedge, 
as  shown  in  Figure  5.2,   it  is  seen  that 

s  =  4)  -  b  -  a  (5.8) 

The  thickness  of  the  freshwater  region,  h,  is  related  to  the  "push  up", 
s,  through 

h  =  6s  -  a  (5.9) 

which  is  a  consequence  of  the  Ghyben-Herzberg  principle. 

The  combination  of  equations  (5.8)  and  (5.9)  results  in 

h  =  6(()  -  6b  -  (1  +  6)a  (5.10) 

But  it  has  already  been  assumed  that  the  relationship  between  h  and  cj) 
in  the  interface  zone  is  given  by  equation  (5.3).  Consequently  the 
values  of  a  and  B  must  be: 

a  =  6   and   3  =  -  6b  (1  +  6)a  (5.11) 

Substituting  these  values  into  equation  (5.4)  gives  the  definition  of 
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the  potential  function  in  the  interface  zone  for  a  confined  coastal 
aquifer, 

$.  =  l&m   -  b  -  (1  +  l/6)a]^  (5.12) 

Now,  in  order  that  the  potential  function  may  be  continuous  throughout 
the  aquifer,  it  is  necessary  to  determine  c  in  equation  (5.5).  Note 
that  at  the  toe  of  the  saltwater  edge 

s  =  (a  +  b)/6  (5.13) 

which  along  with  equation  (5.8)  gives 

({,  =  (a  +  b)  (1  +  1/6)  (5.14) 

Using  this  value  of  ^   in  both  equations  (5.5)  and  (5.12)  and  requiring 
that  they  be  equal  at  the  toe  of  the  wedge  results  in 

c  =  ^  -  Kb(a  +  b)  (1  +  1/6)  (5.15) 

The  definition  of  the  potential   function  in  the  inland  zone  of  a 
confined  coastal   aquifer  then  becomes 

$  =  Kb[<})  +  1/2  b/6-   (a  +  b)    (1   +  1/6)]  (5.16) 

At  the  toe  of  the  saltwater  wedge  this  reduces  to 

^^  =  1/2  Kb^/6  (5.17) 

In  order  to  construct  the  complex  potential  function  for  flow  in 
the  aquifer  the  boundary  condition  at  the  sea  coast  must  be  known. 
From  Figure  5.2  it  is  seen  that  s  =  a/6  at  the  coast.  Substituting  this 
into  equation  (5.12)  gives  the  boundary  condition 
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$  =  0  at  X  =  0  (5.18) 

This  condition  can  be  satisfied  by  locating  an  image  well,  with  the 
same  strength  as  the  injection  well  but  opposite  sign,  on  the  other  side 
of  the  boundary  at  the  same  distance  from  it,  x^^,  as  the  injection  well. 
The  superposition  of  the  velocity  potential  for  the  two  wells  and  the 
uniform  seaward  flow  results  in 

^  =  -xQ    +^}n      I ^=4 A  (5.19) 

((x  +  x^)2+y2) 

When  particular  values  of  ^  are  computed  from  equation  (5.12)  and  subs- 
tituted into  equation  (5.19)  which  is  then  solved  for  the  particular 
equipotential  line,  the  shape  of  the  freshwater  lens  produced  by  the 
injection  well  is  outlined.  Such  a  case  is  shown  in  Figure  5.3. 

From  a  practical  point  of  view  it  is  most  important  to  locate  the 
toe  of  the  saltwater  wedge  around  the  periphery  of  the  lens  because  this 
is  the  outer  boundary  of  the  usable  portion  of  the  freshwater.  The 
equipotential  line  which  delineates  this  outer  boundary  has  the  value 
given  by  equation  (5.17).  Solutions  for  this  line  can  be  made  dimen- 
sionless  by  defining  the  following  dimensionless  parameters: 

(5.20a) 
(5.20b) 
(5.20c) 
(5.20d) 

A  series  of  such  dimensionless  solutions  for  various  conditions  of 
injection  is  shown  in  Figure  5.4.  and  Figure  5.5. 
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Phreatic  coastal  aquifer.  The  Girinskii  potentials  defined  by 
Strack  for  interface  flow  can  also  be  used  in  aquifers  which  have  no 
upper  confining  layer  because  the  thickness  of  the  freshwater  region  is 
still  a  linear  function  of  the  piezometric  head.  The  linear  function 
is  not  the  same  as  for  a  confined  aquifer,  however,  so  the  potentials 
must  be  redefined  to  conform  to  the  aquifer  geometry  shown  in  Figure 
5.6.  In  the  inland  zone  the  specific  discharge,  Q  is  given  by 

The  potential  for  this  region,  then,  should  be  defined  as 

$  =  1/2  k/  +  c  (5.22) 

In  the  interface  zone,  note  that 

s  =  4)  -  b  (5.23) 

and  from  the  Ghyben-Herzberg  principle  the  thickness  of  the  flow  region 
is 

h  =  s(l   +  6)  (5.24) 

which  leads  to 

h  =  (1   +  6)(1)  -   (1   +  6)b  (5.25) 

If  this  is  to  be  in  the  form  of  equation  (5.3)  then  the  constants  a  and 
6  must  be 

a  =  (1  +  6)   and   3  =  -(1  +  6)b  (5.26) 

Substituting  these  values  into  equation  (5.4)  results  in  a  definition 
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for  the  velocity  potential  in  the  interface  zone 

f^.  =  1/2  (1  +  6)  K((})  -  b)^  (5.27) 

At  the  toe  of  the  saltwater  wedge  the  piezometric  head  is 

$^  =  b(l  +  1/6)  (5.28) 

Substituting  this  into  equations  (5.22)  and  (5.27)  and  requiring  that 
be  single  valued  at  this  point  leads  to  the  value  of  the  constant  c  in 
equation  (5.22), 

c  =  -1/2  Kb^  (1  +  1/6)  (5.29) 

The  definition  of  the  velocity  potential  in  the  inland  zone  of  the 
phreatic  aquifer  .then  becomes 

$  =  1/2  K[(j)^  -  (1  +  l/6)b^]  (5.30) 

At  the  toe  of  the  saltwater  wedge  this  reduces  to 

$^  =  1/2  Kb^d  +  &)/S^  (5.31) 

which  is  greater  than  the  corresponding  value  for  a  confined  aquifer  by 
a  factor  of  (1  +  1/6). 

Axi symmetric  Flow  in  a  Leaky  Aquifer 
Consider  a  homogeneous  saline  aquifer  with  hydraulic  conductivity, 
K,  of  uniform  thickness,  b,  and  unlimited  areal  extent.  It  is  confined 
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at  the  bottom  by  a  horizontal  impervious  layer  and  at  the  top  by  a 
horizontal  aquitard  of  thickness,  b',  and  hydraulic  conductivity,  K'. 
A  fully  penetrating  well  injects  freshwater  at  a  constant  rate,  Q, 
forming  a  freshwater  lens,  which  remains  centered  on  the  well  since 
there  is  no  pre-existing  flow  in  the  surrounding  salt  water.  Injection 
continues  until  a  steady  state  is  eventually  reached  in  which  the  injec- 
tion rate  is  matched  by  the  rate  of  leakage  through  the  aquitard.  The 
mixing  zone  at  the  boundary  of  the  lens  then  becomes  relatively  thin 
and  can  be  approximated  by  a  sharp  interface.  This  situation  is 
depicted  in  Figure  5.7. 

The  flow  region  is  divided  into  two  zones,  each  of  which  has  its 
own  governing  equation.  In  zone  1,  which  is  occupied  entirely  by 
freshwater,  the  customary  analysis  for  well  flow  in  leaky  aquifers  can 
be  used.  This  includes  the  assumption  of  90  degree  streamline  refrac- 
tion at  the  aquitard  so  that  flow  in  the  aquifer  is  horizontal  while  in 
the  aquitard  it  is  vertical.  The  governing  equation  is  (see  De  Weist, 
1965,  p.  264) 

§^1^-  (1)^  =  0  (5.32, 


where 


B^  =  Kbb'/K' 


This  is  a  modified  Bessel  equation  of  order  zero,  for  which  the 
general  solution  is  given  by 

s  =  C^jQ(r/B)  +  C2i^Q(r/B)  (5.33) 

where  Jq  and  K^   are  the  modified  zeroth  order  Bessel  functions  of  the 
first  and  second  kinds,  respectively.  The  boundary  condition  at  the 
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outer  edge  of  zone  1  is  complicated  by  the  fact  the  position  of  the 
boundary  is  unknown.  Since  the  Ghyben-Herzberg  relationship  is  being 
used,  the  value  of  s  at  the  toe  of  the  interface  is  known,  though,  and 
this  will  be  used  to  determine  the  boundary  location.  The  boundary 
conditions  for  zone  1,  then,  are 


^)  ri\     ^'P   =  -   Q/(2TrKb)  (5.34a) 


2)  s  =  ^  at  r  =  r^  (5.34b) 


Differentiating  equation  (5.33)  and  applying  the  first  boundary  condi- 
tion gives 

Ci  ^^i'"o  [""/B  J/r/B)]  -  C^  ^'^i\   [r/B  /f^(r/B)]  =  -Q/(2TTKb)    (5.35) 

in  which  i^  and  K-^  are  modified  Bessel  functions  of  order  1.  The  first 
limit  in  equation  (5.35)  is  zero  and  the  second  is  one,  which  leads  to  a 
value  for  the  constant,  Co. 

Substituting  this  and  the  second  of  equations  (5.34)  into  equation  (5.33) 
gives  the  other  required  constant 

S  =  t^-  2^^0(^t/B)]/^0(^t/B)  (5.37) 

The  solution  for  s  in  zone  1  is  then 

rb+a    0  ^r)ir/^)  n 

in  which  the  value  of  r.  is  as  yet  undetermined. 
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The  flow  in  zone  2  is  complicated  by  the  presence  of  the  interface 
so  that  a  new  equation  must  be  derived  for  it.  In  this,  the  Dupuit 
assumption  will  be  used  so  that,  as  before,  only  head  loss  due  to 
horizontal  flow  in  the  aquifer  will  be  considered. 

Applying  the  principle  of  mass  conservation  to  the  differential 
wedge  of  aquifer  shown  in  figure  5.8  gives 

uhrdo  -  (u+du)(h+dh)(r+dr)de  -  w(s)rdrde  =  0  (5.39) 

expanding  this,  cancelling  terms  where  possible,  and  dividing  by  drde 
results  in 

^•^  dr  "^  ^"^  ^  "^  "^  3?  ^"^  "^  ^h  +  dhu  +  dudh  +  rw(s)  =  0       (5.40) 
Taking  the  limit  as  the  differentials  go  to  zero,  this  becomes 


Substituting  Darcy's  law 


"  =  -  K^  (5.42) 


into  equation  (5.41)  and  dividing  by  Khr  gives 

d^  .  1  dh  ds  ,  1  ds   w(s)  _  _  ,    , 

^^2   h  ^  dr   r  dr  "  UT  ~   "  (5-43) 

But  the  Gyben-Herzberg  principle  gives  the  following  relationship 
between  h  and  s 

h  =  (6s-a)  (5.44) 
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Substituting  this  into  equation  (5.43)  gives 


2 
d  s  .   6   ,dsx2  ,  1  ds   w(s) 
'TJ^W^'^d^^     ^rd^--W-=^  (5.45) 


Noting  that 


w{s)  =  K's/b'  (5.46) 

the  governing  equation  for  flow  in  zone  2  is  obtained 

d^s    6    ,ds^^  ,  1  ds   K-s     ^  ,,  _, 

^^2^l6i3i7W^  ""Fd?-  (6s-a)b'K  =  °  (5.47) 

Equation  (5.47)  is  non-linear  and  is  not  one  of  the  few  types  of 
non-linear  equations  for  which  analytical  methods  of  solution  are 
available.  There  are  several  routine  numerical  methods  which  can  handle 
such  an  equation,  but  first  it  is  necessary  to  find  the  correct  boundary 
conditions  by  which  the  flow  in  zone  2  can  be  coupled  to  that  in  zone  1. 
Those  conditions  are  partially  provided  by  the  requirement  that  the 
piezometric  head,  s,  and  its  slope,  ^  ,  must  be  continuous  at  the  radius, 
r^,  where  the  two  zones  join.  Now,  the  value  of  s  at  that  point  has 
already  been  given,  and  its  slope  can  be  found  by  differentiating  equa- 
tion (5.38).  These  requirements  can  be  stated  as: 

1)  ^(^t)  =^  (5.48a) 


2)  ^ 

'     dr 


rb+a    0  •^l('"t/^^    n 

=  ^~  -  2^^0(V^)^  BjJr,/B)  -  2^^l(VB)  (5-48b) 
t  '^     ^ 


These  are  only  partial  boundary  conditions,  however,  because  they  con- 
tain the  unknown,  r^.  One  final  boundary  condition  is  needed  which  can 
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be  described  as  the  requirement  that  as  r  approaches  infinity,  s  and  all 
of  its  derivatives  should  approach  zero.  The  truth  of  this  is  not 
immediately  apparent,  but  it  will  be  demonstrated  after  some  of  the 
integral  curves  to  the  set  of  coupled  equations  (5.32)  and  (5.47)  have 
been  obtained. 

The  integral  curves  can  be  obtained  by  choosing  values  for  the 
variable,  r^.  When  this  is  done,  the  solution  for  zone  1  is  immediately 
obtained  from  equation  (5.38).  The  boundary  conditions  (5.48)  then 
provide  the  necessary  starting  point  for  a  numerical  treatment  of  equa- 
tion (5.47).  The  integral  curves  shown  in  Figure  5.  9  were  obtained  by 
the  Adams-Bashforth  type  predictor-corrector  method  using  the  Runge- 
Kutta  method  as  a  starter  for  the  first  four  distance  steps.  Figure 
5.9  reveals  two  different  kinds  of  integral  curves  representing  solu- 
tions for  equation  (5.47)  in  terms  of  s  vs^.  r  for  different  values  of 
r^.  In  the  upper  curves,  s  first  starts  to  decrease  with  increasing  r, 
then  a  minimum  value  is  reached  and  s  starts  to  increase.  This  is 
obviously  not  a  physically  possible  solution  to  the  problem.  The  lower 
curves,  on  the  other  hand,  continue  to  approach  zero  until  they  reach 
it.  This  is  not  a  physically  possible  steady  state  solution  either 
because  as  s  approaches  zero  along  this  line,  its  slope  becomes  increas- 
ingly negative  which  means  that  the  flow  velocity  is  increasing.  If 
this  were  true,  it  would  mean  that  the  flow  field  engendered  by  the 
injection  well  has  a  finite  boundary  and  that  the  flow  rate  at  that 
boundary  is  ^jery   high.  It  is  instructive  here  to  consider  the  equation 
for  the  curvature  of  the  piezometric  head  line  in  the  special  case  of 
a  =  0. 

^-     ^'  ■>  fdsi^   1  ds  ,,  .-, 


90 


«3- 


in 


ex: 
o 


LU  to 


< 


en 
to 

UJ 

o 


(Sd3i3W)    s 


91 


This  shows  that  as  s  approaches  zero  with  non-zero  slope,  the  curvature 
becomes  negatively  infinite.  This  is  what  is  happening  in  the  lower  set 
of  integral  curves.  However,  if  both  s  and  ^  approach  zero  in  the 
right  way,  the  curvature  ^will  approach  zero  also.  This  happens 
somewhere  in  the  area  between  the  upper  and  lower  integral  curves,  and 
it  represents  the  only  physically  possible  solution.  While  it  is  not 
possible  to  exhibit  this  solution  precisely,  it  is  possible  to  obtain 
the  upper  and  lower  integral  curves  as  close  to  it  as  one  may  desire. 

Figure  5.10  shows  solutions  for  the  toe  of  the  interface,  r. , 
which  were  claculated  by  this  method  for  a  variety  of  flow  conditions. 
They  are  presented  in  the  form  of  dimensionless  variables  6Q/(2TrKb^)  vs^. 
r^/B  for  several  values  of  the  density  ratio,  a/b,  on  each  figure.  Of 
course,  knowledge  of  the  value  of  r^  does  not  constitute  a  complete 
solution  to  the  problem.  One  would  also  like  to  know  the  shape  of  the 
interface  and  consequently  of  the  entire  lens.  Examples  of  such  solu- 
tions are  shown  in  Figure  5.11.  It  should  be  noted,  however,  that  from 
a  practical  standpoint,  the  position  of  the  toe  of  the  interface  is  most 
important,  since  the  freshwater  located  in  that  part  of  the  lens 
beyond,  r^,  cannot  be  recovered  unmixed  with  salt  water.  In  fact,  it 
will  not  be  possible  even  to  recover  all  of  the  water  inside  r.  because 
of  dispersion  and  further  tilting  of  the  interface  when  the  lens  departs 
from  its  equilibrium  position. 
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CHAPTER  6 
A  FINITE  ELEMENT  MODEL  FOR  THE  STEADY  FLOW  CASE 

General 


The  solutions  to  steady  state  problems  given  in  the  preceeding 
chapter  depended  on  the  assumption  of  horizontal  streamlines  and  static 
conditions  in  the  saline  part  of  the  aquifer.  These  two  simplifications 
were  required  so  that  the  Ghyben-Herzberg  principle  could  be  used  for 
location  of  the  position  of  the  interface.  If  the  horizontal  streamline 
restriction  is  lifted,  a  three-dimentional  boundary  value  problem  must 
be  solved  and  the  interface  must  be  located  by  a  trial  and  error 
process.  If  flow  is  permitted  in  the  saline  portion  of  the  aquifer  as 
well,  then  a  pair  of  coupled  three-dimensional  boundary  value  problems 
must  be  solved.  In  this  chapter  the  method  of  finite  elements  will  be 
used  to  obtain  numerical  solutions  to  such  problems. 

Description  of  the  Galerkin  Finite  Element  Scheme 
The  Fundamental  Concept 

The  boundary  value  problem  within  a  steady  freshwater  lens  can  be 
posed  in  terms  of  the  piezometric  head  cf)  =  (})(x,  y,  z),  a  dependent 
variable  which  is  a  continuous  function  of  position  inside  the  bounda- 
ries of  the  lens.  In  the  finite  element  method  this  continuous  function 
is  approximated  by  a  set  of  piecewise  continuous  polynomial  functions, 
each  of  which  is  defined  on  a  discrete  subregion,  or  element,  of  the 
flow  domain.  In  order  to  do  this,  the  flow  region  must  first  be 
divided  up  into  elements  with  a  finite  number  of  nodal  points  denoting 
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the  boundaries  of  the  elements.  The  values  of  the  continuous  quantity, 
(j),  at  all  of  these  nodal  points  form  a  set  of  unknown  variables  for  which 
the  boundary  value  problem  is  to  be  solved.  The  piecewise  continuous 
polynomial  functions  are  defined  over  the  elements  in  terms  of  these 
unknown  nodal  values  of  ^.     For  this  reason,  they  are  called  interpola- 
tion functions. 

Once  the  set  of  interpolating  functions  is  obtained  the  values  of 
the  nodal  unknowns  must  be  determined  in  such  a  way  that  the  interpola- 
tion polynomials  will  provide  the  closest  possible  approximation  to  the 
governing  differential  equation.  There  are  several  ways  of  doing  this. 
The  method  used  here  is  to  minimize  the  error  of  the  approximating 
functions  when  integrated  over  the  flow  domain  by  the  Galerkin  weighted 
residual  technique. 

The  Element  and  its  Shape  Functions 

The  type  of  element  to  be  used  is  illustrated  in  Figure  6.1.  It 
is  a  three-dimensional  element  with  20  nodes.  The  presence  of  the  mid- 
side  nodes  on  each  side  permits  the  dependent  variable,  4),  to  be 
approximated  by  an  interpolation  function  that  is  quadratic  in  all  three 
coordinate  directions.  The  form  of  this  interpolation  polynomial  is 

(t>  =  a^  +  a^  X     +  a^  y     +0^2  +  a^  xy  +  05  xz  +  a^  yz 
+  Ug  xyz  +  ag  x2  +  a^Q  y^  +  a-^  z^  +  a^^  x^y     +  a^^  x^z 
+   a^4  xy2  +  a^^  y2z  +  a^^  xz^  +  a^^  yz2  +  a-,g  x^yz 
+  a^g   xy2z  +  a^Q  xyz^  (6.1) 

where  (p   is  the  approximation  to  0  and  the  20  coefficients  a,  .  .  .  a^^ 
must  be  defined  in  terms  of  the  element  shape  and  the  unknown  values  of 
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(J)  at  the  20  nodes,  O^  .  .  .  ^^q-  The  definition  of  these  coefficients 
must  be  done  in  such  a  way  that  the  values  of  (^  given  by  equation  (6.1) 
will  satisfy 

*  (x,  y,  z)  =  <D.   when   (x,  y,  z)  =  (X.,  Y • ,  Z . )     (6.2) 

where  X^. ,  Y^.  and  Z^  are  the  coordinates  of  the  i^*^  node.  It  is  possible 
to  apply  the  condition  of  equation  (6.2)  to  equation  (6.1)  for  all  20 
nodes  and  solve  the  resulting  20  equations  for  a,  .  .  .  a^Q.  Substitu- 
ting the  results  of  this  into  equation  (6.1)  it  can  be  rearranged  into 
the  following  form: 

^   20 

<i>  =    lu    N.  $.  (6.3) 

i=l    ^ 

where  the  20  coefficients  N.  •  •  •  N20  are  called  shape  functions 
because  they  depend  only  on  the  geometry  of  the  element.  Since  equation 
(6.3)  must  also  satisfy  the  requirement  of  equation  (6.2)  it  is  obvious 
that  at  the  i  node  N^.  =  1  and  all  the  other  N's  are  zero.  The  values 
of  the  shape  functions  between  nodes  must  be  such  that  cfi  will  be  conti- 
nuous over  the  element.  The  distributions  of  two  quadratic  shape  func- 
tions over  a  two-dimensional  element  are  illustrated  in  Figure  6.2. 

The  element  and  its  shape  functions  can  be  considered  to  adequately 
represent  the  flow  domain  and  the  dependent  variable,  0,  if  they  result 
in  a  description  which  converges  to  the  true  distribution  of  <t>   as  the 
flow  domain  is  subdivided  into  finer  and  finer  elements.  It  can  be 
shown  (Zienkiewicz,  1971)  that  the  solution  will  be  convergent  if  two 
conditions  are  satisfied.  First,  the  element  shape  functions  must  be 
capable  of  correctly  describing  a  constant  distribution  of  the  dependent 
variable  throughout  the  element.  If  the  dependent  variable  were  constant 
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a.  Distribution  of  shape  function  N, 


b.  Distribution  of  shape  function  Uj 


FIGURE  6.2.  DISTRIBUTION  OF  TWO  QUADRATIC  SHAPE  FUNCTIONS  OVER 
A  TWO-DIMENSIONAL  QUADRILATERAL  ELEMENT 
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then  all   twenty  values  of  <J>p  (t^,   .    .    .   <t>^Q  would  be  equal  and  ^  would 
be  correctly  described  if 

20 

E     ^i    (x,  y,   z)   -  1  (6.4) 

1=1 

for  all  values  of  x,  y,  and  z  in  the  element.  Second,  the  predicted 
values  of  the  dependent  variable,  <t>,   must  be  continuous  across  the 
boundaries  between  adjacent  elements.  In  order  that  this  may  be  assured, 
it  is  necessary  that  the  shape  function  values  for  the  adjacent  elements 
must  be  equal  along  their  common  boundaries.  The  20  node  element  shown 
in  Figure  6.1  is  one  of  a  family  called  the  Serendipity  family  of 
elements  because  20  shape  functions  can  be  found  for  it  by  inspection 
which  satisfy  all  of  the  above  requirements. 

The  20  shape  functions  are  defined  in  terms  of  the  local  coordinate 
system  shown  in  Figure  6.3.  The  adoption  of  a  local  coordinate  system 
is  necessary  for  the  quadrilateral  family  of  elements  to  insure  inter- 
element  continuity  when  the  element  boundaries  are  not  parallel  to  the 
global  X,  y,  z  coordinates  (Segerlind,  1976).  The  shape  functions  are 
given  in  Table  6.1 . 

Finite  Element  Formulation  of  the  Governing  Equation 

The  governing  equation  for  steady  groundwater  flow  in  a  homogeneous 
leaky  aquifer  is 

M*)  =  K  ^^  4  ^  ^y  4  'hz^    +  Q  -  ^  =0     (6.5) 
where  Q  is  the  strength  of  a  fluid  source  at  the  surface  of  the  flow 

2  9 

region,  g  =  b'/K'  is  a  modified  form  of  the  leakage  coefficient,  B^, 
and  (pQ   is  the  piezometric  head  in  the  phreatic  aquifer  above  the 
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TABLE  6.1 
SHAPE  FUNCTIONS  FOR  THE  20  NODES  ELEMENT 

N^  =  1/8  (1  -  d  (1  -  n)  (1  -  d  (-  ?  -  n  -  C  -  2) 
N,  =  1/8  (1  +  c)  (1  -  n)  (1  -  ?)  (?  -  ri  -  2;  -  2) 


^3 

^'4 


^5 


1/8  (1  +  0  (1  -  n)  (1  -  d  (C  +  n  -  c  -  2) 

1/8  (1  -  ^)  (1  +  n)  (1  -  d  (-  5  +  n  -  c  -  2) 

1/8  (1  -  ^)  (1  -  n)  (1  +  d  (-  C  -  n  +  C  -  2) 

1/8  (1  +  C)  (1  -  n)  (1  +  c)  (C  -  n  +  ^  -  2) 


N.  =  1/8  (1  +  0  (1  +  n)  (1  +  C)  (C  +  n  +  C  -  2) 


N„  = 


'10 


U 


N 


N 


^11 
^12 
13 
14 
15 
16 
^17 
^18 
^19 
^20 


1/8  (1 
1/4  (1 
1/4  (1 
1/4  (1 
1/4  (1 
1/4  (1 
1/4  (1 
1/4  (1 
1/4  (1 
1/4  (1 
1/4  (1 
1/4  (1 
1/4  (1 


0  (1  +  n)  (1  +  ?)  (-  C  +  n  +  C  -  2) 


C)  (1  -  n 

n^)  (1  +  C 

K^)  (1  +  n 

n^)  (1  -  C 


.2)(1 


-  c')  (1  +  ? 

-  C^)  (1  +  ? 

-  C^)  (1  -  <^ 

-  C^)  (1  -  n 


-  n^)  (1  +  K 


-   r)  (1  +  n 


n")  (1  -  ? 


(1 
(1 
(1 
(1 


(1  -  n) 

(1  -  n) 

(1  +n) 

(1  +  n) 

(1  +  rj 

(1  +  ?) 

(1  +  0 

(  +  0 
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leaky  layer.  As  described  in  a  previous  section,  the  finite  element 
method  is  to  approximate  the  variable  ())  by  the  piecewise  continuous 
function  (|).  If  this  could  be  done  perfectly  then  equation  (6.5)  could 
be  written  in  terms  of  (j)  as 

L  U)   =  0  (6.6) 

for  every   point  in  the  flow  domain.  In  reality,  the  best  that  can  be 
done  is  to  minimize  the  operation  L   (cj))  in  an  average  sense  over  e\/ery 
element  by  the  Galerkin  weighted  residual  method,  so  that 


/ 


N.  L  ((J>)  dV  =0;  for  j  =  1,  2,  .  .  .  20        (6.7) 
in  which  the  weighting  functions  N.  are  the  same  as  the  element  shape 

J 

functions.  In  expanded  form  this  is 
/  kx  "i  ^  -vv  N,  4  ^  K„  N,  4  .  N,Q  -  N,  '*:^'ldV=0  ,6.8, 


'V 


9x     ^^  ^  9y"    ^^  ^  9z^    -J     J   3' 


It  should  be  recalled  that  the  element  shape  functions  listed  in 
Table  6.1  were  required  to  satisfy  the  condition  that  only  the  value  of 
(})  need  be  continuous  between  elements.  Since  the  continuity  of  the 
first  derivative  of  <t)   is  not  guaranteed  the  integral  in  equation  (6.8), 
which  contains  second  derivatives  of  0,  may  not  necessarily  exist.  For 
this  reason  the  order  of  the  derivaties  in  equation  (6.8)  must  be 
reduced  by  one.  This  is  done  through  the  use  of  the  divergence  theorem 
and  the  product  rule  for  derivatives  which  states  that 

9     9(})    9N.  9(J)      9^^ 

—  (N.  — )  =  -J-— +  N.  -^  (6.9) 

3x   "J  ax     9x  3x    -^   9x 

Rearranging  this  and  substituting  into  equation  (6.8)   results  in 
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K....  ^  (N,  U)   +  K 


'  (N,  |i)  .  K     A  (,.  11)1  dV 


y    t  XX   9x   ^"j    3x'        >y   8y   ^'^j    Sy^   "  ^z   8?  ^'^j    8^^^ 


i  h ' 


{(f)-(J)Q)   N./6 


J       ) 


dV 


/    [k       9'^i   9*  3N.    9(f)  3N.   9 J 


9x    9x        yy    9y    9y       "^zz    9z    9z^ 


dV 


(6.10) 


Noting  that,  by  Darcy's  law, 


\x  8x  ^"^j   9x^ 


^  (^  ^x) 


(6.11) 


the  first  integral  in  equation  (6.10)  can  be  changed  by  the  divergence 
theorem  into 


X 


N,   q    •    n  dA 


when  n  is  the  outward  unit  vector  normal   to  the  boundary,  S. 

Now  consider  the  terms  in  the  second  integral  of  equation  (6.10). 
The  first  term,  containing  Q,  refers  to  the  addition  of  fluid  at  the 
bounding  surface  of  the  flow  domain.     The  second  refers  to  fluid  lost 
from  the  aquifer  by  leakage  through  the  boundary  in  contact  with  the 
leaky  layer.     The  second  integral,   then,  can  be  changed  to  a  surface 
integral   because  the  integral   is  zero  everywhere  except  on  the  surface 
of  the  flow  region.     Equation  (6.10)   then  becomes 


/ 


I         9N.    9(J, 
k,...  —i—  +  K 


9N.   d<p 


+  K. 


9N.   U] 


xx    9x    9x       >y   9y    9y       ^zz    9z    9z 


■/ 


Jj  q   •   n  +  NjQ  -   (^-(^q)   H./B^ 


dV 


dA 


(6.12) 
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It  was  shown  in  the  previous  section  that  the  approximated  piezo- 
metric  head,  (j),  is  defined  over  each  element  in  terms  of  the  element 
shape  functions  and  the  nodal  unknowns  by  equation  (6.3).  Applying  this 
definition  to  equation  (6.12)  gives  the  finite  element  formulation  of 
the  governing  equation  for  an  element 


■./ 


K 


9N.    3N. 


+  K 


9N.    9N. 


+  K. 


9N.    9N.1 


XX    9x     9x         yy    9y     9y         zz    9z     9z 


dV 


■i  J   ^•^■/^^  "^^  ^    J    [Nj  q   •    n  +  NjQ  +  Nj^B^j   dA       (6.13) 


for 

i   =  1,  2, 
J  =  1,  2, 

In  matrix  form  this  can  be  written  as 


20 
20 


[k]  {<!'}=  {f} 


(6.14) 


where  [k]  is  a  20  x  20  matrix  known  as  the  element  stiffness  matrix,  {f}, 
is  a  column  matrix  of  forcing  functions  given  by  the  right  hand  side  of 
equation  (6.13),  and  {$}  is  the  matrix  of  nodal  unknowns,  ^■.,   ^2^,   .    .  . 

'^20- 

In  accordance  with  the  fundamental  concept  of  the  finite  element 

method,  a  similar  equation  can  be  written  for  the  entire  flow  domain 

by  a  proper  summation  of  matrix  equations  for  the  individual  elements. 

This  results  in 


[K]  W  =   {F} 


(6.15) 
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in  which  [K]  is  the  global  stiffness  matrix,  {F}  is  the  global  matrix 
of  forcing  functions,  and  W   is  the  matrix  of  unknown  values  of  piezo- 
metric  head  at  all  nodes.  These  matrices  are  obtained  from  the  element 
matrices  by  standard  methods  of  assembly  which  can  be  found  in  any  finite 
element  textbook. 

Development  of  the  Element  Matrices 

The  terms  of  the  element  stiffness  matrix  are  given  by  the  inte- 
grals on  the  left  hand  side  of  equation  (6.13)  in  which  the  element 
shape  functions  are  prominent.  The  evaluation  of  these  integrals  is 
complicated  by  the  fact  that  the  element  shape  functions  are  defined  in 
the  local  element  coordinate  system  while  the  governing  equation  refers 
to  the  global  x,  y,  z  system.  The  transformation  of  coordinate  systems 
is  accomplished  by  use  of  the  Jacobian  matrix,  [J],  defined  as 


[J]  = 


[ 


3x 

95 

9y 

9? 

3z' 

95 

8x 

9n 

9ti 

9z 

9n 

3x 
9? 

9y 

9c 

9z 
9cJ 

(6.16) 


The  terms  in  this  matrix  depend  on  the  geometry  of  the  element  in  such 
a  way  that 

20  dU, 


'^-  Z  w  ^• 


95 


(6.17) 


i=l 


so  that  the  Jacobian  matrix  for  an  element  becomes 
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[J]  = 


20 

8N. 

20 

9N. 

20 

9N. 

L 

8C 

^• 

E 

1 

9C 

^• 

E 

1 

9e; 

1=1 

i=l 

i=l 

20 

9N. 

20 

9N. 

20 

9N. 

E 

9n 

^• 

E 

1 

9n 

^• 

E 

1 

9n 

i  =  l 

i=l 

i=l 

20 

9N. 

20 

9N. 

20 

9N. 

E 

1 

9C 

^• 

E 

1 

^• 

E 

1 

9C 

i  =  l 

i=l 

i=l 

(6.18) 


It  should  be  remarked  here  that  the  representation  of  the  relation- 
ship between  the  two  coordinate  systems  given  by  equation  (6.17)  is  not 
the  only  possible  one.  If  the  boundaries  of  the  element  were  restricted 
to  planar  surfaces  only,  there  would  be  no  need  to  use  the  coordinates 
of  all  twenty  nodes  in  the  description  of  element  shape.  Only  the 
corner  nodes  would  be  essential  because  they  give  enough  information  for 
the  location  of  planar  boundaries.  In  such  a  case  the  element  would  be 
called  subparametric  because  the  degree  of  the  equations  defining  the 
element  shape  is  less  than  the  degree  of  the  interpolation  polynomials. 
The  relationship  given  by  equation  (6.17)  permits  the  element  boundaries 
to  be  curved  surfaces  described  by  quadratic  equations.  Since  these 
equations  are  of  the  same  degree  as  the  interpolation  polynomials  this 
is  called  an  isoparametric  element. 

The  Jacobian  matrix  is  used  in  two  ways  in  the  evaluation  of  the 
first  integral  of  equation  (6.13).  First,  the  integrand  calls  for  the 
evaluation  of  partial  derivatives  of  the  shape  functions  with  respect 
to  the  X,  y,  z  coordinate  system.  This  can  be  accomplished  by  noting 
from  the  chain  rule  for  derivatives  that 
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3N-   9N.  8x   9N.  dy       8N.  8z 
—  -       +  -~-  ^^  +   ^ 


3C 


9x  3C    9y  3^    9z  9? 


(6.19) 


and  that  9x/9C,  9y/9C,  and  9z/9?  are  terms  of  the  Jacobian  matrix. 
Consequently,  we  have 


9N^. 


95 

9N. 
"9^ 

9N. 

"9r 


=  [J] 


9N^. 
Ix" 

9N. 
■97 

9N. 


(6.20) 


and 


9N. 
"97 


9N. 


9N. 


[J]' 


9N. 

9N. 

"9n" 


9N. 

"9r 


(6.21) 


The  expression  for  the  partial  derivatives  with  respect  to  the  local 
coordinate  system  are  given  in  Table  6.2. 

The  second  use  of  the  Jacobian  concerns  the  evaluation  of  the 
volume  integral  itself.  The  integral  will  be  evaluated  in  terms  of  the 
element's  local  coordinate  system,  in  which  the  values  of  each  coordinate 
vary  between  +  1  and  -  1 .  In  order  to  give  the  integral  the  correct 
volume  the  change  of  scale  included  in  the  coordinate  transformation 
must  be  accounted  for.  This  is  done  by  multiplying  the  integrand  by 
the  absolute  value  of  the  determinant  of  the  Jacobian  matrix,  Idet  J|. 
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TABLE  6.2 
PARTIAL  DERIVATIVES  OF  THE  SHAPE  FUNCTIONS 


9N^/9?  =  1/8  (1  -  n)  (1  -  ?)  (2?  +  n  +  ?  +  1 

m^/^K  =  1/8  (1  -  n)  (1  -  0  (2^  -  n  -  c  -  1 

dU^/di  =  1/8  (1  +  n)  (1  -  c)  (2^  +  n  -  c  -  1 

dN^/dK  =  1/8   (1   +  n)  (1    -   c)  (2^  -  n  +  ?  +  1 

SNg/ac  =  1/8  (1  -  n)  (1  +  c)  (2c  +  n  -  c  +  1 

sNg/s^  =  1/8  (1  -  n)  (1  +  d  (2C  -  n  +  c  -  1 

aN^/ac  =  1/8  (1  +  n)  (1  +  ?)  (2c  +  n  +  ?  -  1 

3Ng/3?  =  1/8  (1  +  n)   (1  +  d   (2^  -  n  -  C  +  1 

dHg/dK  =  -  1/2C  (1  -  n)   (1  -  c) 
9N^0/9C  =  1/4   (1    -  n^)    (1-   C) 
3N^^/3C  =  -   1/2^   (1   +  S)    (1    -   C) 
3N^2/9^  =  -   9N^q/3? 
3N^3/35  =  -   1/4   (1    -   C^)    (1    -  n) 
3N^4/3^  =   3N^3/3^ 
3N^5/3?  =   1/4   (1    -   ^2)    (1    +  n) 

3N^g/35  =  aN^g/ac 

3N^7/3C  =  -   1/2C   (1    -  n)    (1   +  n) 
3Ni8/9^  =  1/4   (1    -  n^)    (1   +  d 
3N^g/3C  =  -   1/2C   (1    +  n)    (1    +  C) 
3N2q/3C  =  -   3N^8/9^ 
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TABLE  6.2  -  continued 


9N^/9ri 


8N2/8n 


9N3/9n 


SN^/Sn 


=  1/8   (1  -   U    (1   -   C)    (2n  +  ?  +  C  +  1 

=  1/8  (1  +  o   (1  -  0   (2n  -  C  +  ?  +  1 

=  1/8  (1  +  d   (1  -  U   (2n  +  ?  -  ?  -  1 

=  1/8  (1  -  a   (1  -  C)   (2n  -  ^  -  c  -  1 


SNg/an 


=  1/8  (1  -  0  (1  +  c)  (2n  +  ?  -  ?  +  1 


9Ng/9n 

3N7/9n 

9Ng/3n  = 

9N9/9T1  = 

9N^Q/9n  = 

9N^^/9n  = 

9N^2/9n  = 

9N^3/9C  = 

9N^4/9n  = 

9^1 5/ 9n  = 

9N^g/9n  = 

9N^7/9n  = 

9N^g/9Ti  = 

9N^g/9n  = 


9N2Q/9n  = 


=  1/8  (1  +  a   (1  +  d   (2n  -  ?  -  S  +  1 


=  1/8  (1  +  U   (1  +  C)   (2n  +  ?  +  c  -  1 
=  1/8  (1   -  ?)   (1   +  c)   {2n  -  e  +  C  -  1 


1/4  (1  -  t)  (1  -  K) 
l/2n  (1  +?)(!-  s) 
9Ng/9n 

l/2n  (1  -  U  (1  -  c) 

1/4  (1  -  c^)  (1  -  c) 

1/4   (1    -   C^)    (1   +  a 

9N^4/9n 

9ri^3/9n 

1/4   (1   -   C^)    (1   +  c) 

l/2n  (1  +  a  (1  +  -;) 

9N^7/9n 

l/2n  (1  -  C)   (1  +  C) 


no 


TABLE  6.2  -  continued 


3N^/9^  =  1/8 
9N2/9C  =  1/8 
dU^/dt;  =  1/8 
8N^/9?  =  1/8 
dH^/dr,  =  1/8 
9Ng/9C  =  1/8 
9N7/9C  =  1/8 
9Ng/9?  =  1/8 


1  -  C)    (1   -  n)  (2?  +  e  +  n  +  1 

1  +  0  (1  -  n)  (2c  -  ^  +  n  +  1 

1  +  0  (1  +  n)  (2c  -  C  -  ri  +  1 

1  -  ^)  (1  +  n)  (2c  +  ^  -  n  +  1 

1  -  ?)   (1  -  n)  (2c  -  ?  -  n  -  1 

1  +  0  (1  -  n)  (2c  +  C  -  n  -  1 

1  +  d  (1  +  n)  (2c  +  ?  +  n  -  1 

1  -  0  (1  +  n)  (2c  -  C  +  n  -  1 


9  Ng/9c  =  -  1/4  (1  -  n"^)  (1  -  d 

9N^q/9c  =  -   1/4   (1   -  n^)    (1   +  0 

9N^^/9C  =  -   1/4   (1   -  K^)    (1   +  n) 

9N^2/9?  =  -   V4   (1    -  n^}    (1    -   K) 

9N^3/9C  =  -   1/2   (1    -  d    (1    -  n) 

9N^4/9C  =  -   1/2   (1   +  d    (1    -  n) 

9N^5/9C  =  -   1/2   (1   +  0    (1    +  n) 

3N^g/9C  =  -   1/2   (1    -  a    (1    +  n) 
9N^7/9C  =  -    9N9/9C 
9N^g/9C  =  -    9N^q/9c 
9N^g/9C  =  -   9N^^/9C 
9N2q/9c  =  -    9M^2/^^ 


m 


For  those  elements  that  do  not  have  leaky  boundaries  the  terms  of  the 
element  stiffness  matrix,  k..,  are  given  by 


^ij 


1  .1  .1 


III 


9N.  9N. 


9N.  9N. 


yy    9y  9y 


+  k 


9N.  9N.1 


ZZ  9Z   9Z  ^ 


|det  J|  ds  dn  d^ 


(6.22) 


for 


1=1,2, 
J  =  1,  2, 


20 
20 


where  the  partial  derivatives  are  obtained  from  equation  (6.21)  and 
Table  6.2. 

Since  the  functions  in  the  integrand  of  equation  (6.22)  depend  on 
the  geometry  of  the  element,  the  integral  must  be  evaluated  numerically. 
Since  the  shape  functions  appearing  in  the  integrand  are  quadratic  in 
?,  n,  and  c  the  integral  can  be  evaluated  exactly  by  Gaussian  quadrature, 
using  the  eight  sample  points  shown  in  Figure  6.4  with  a  weighting  factor 
of  1.0  for  each  point.  Equation  (6.22)  then  becomes 


8 


n=l 


where 


(6.23) 


g  =  the  integrand  of  equation  (6.22) 
^n'  "^n'  ^n  ^  ^^^   coordinates  of  the  n   sample  point 
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If  the  element  contains  a  leaky  boundary  the  terms  of  the  element 
stiffness  matrix  will  also  have  contributions  from  the  second  integral  in 
equation  (6.13).  This  is  a  surface  integral  and,  since  leakage  only 
occurs  through  surfaces  in  contact  with  the  leaky  layer,  it  is  only 
integrated  over  the  top  side  of  the  element.  If  the  nodes  of  all  ele- 
ments are  numbered  consistantly  as  shown  in  Figure  6.1  then  the  only 
shape  functions  contributing  to  the  surface  integral  are  those  associated 
with  nodes  5,  6,  7,  8,  17,  18,  19,  and  20.  The  others  are  all  zero  on 
this  surface.  The  integral  then  becomes 

/  N.Nj/3^  ^^  ^     J      J    "^i^j/^^Net  J I  dn  d^        (6.24) 

for 

i  =  5,  6  ...  8,  17,  18  ...  20 
j  =  5,  6  ...  8,  17,  18  ...  20 

Gaussian  quatrature  for  this  surface  integral  requires  the  nine  sample 

points  shown  in  Figure  6.5  and  their  respective  weighting  factors.  The 

surface  integral  is  evaluated  numerically  by 

J    1  9 

J        J      N.Nj./a^ldet  J|  dn  dC  =  E  W^  glc^,  n^)       (6.25) 


where 

g  =  the  integrand  of  the  integral  on  the  left  hand  side 
C  ,  n_  =  the  coordinates  of  the  n   sample  point 

W  =  the  weighting  factor  for  the  n   sample  point 

The  element  stiffness  matrix  for  a  leaky  element  is  obtained  by  adding 
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8 


Weight 

Sample 

Factor 

Point 

n 

.    "n 

1 

-a 

-a 

25/81 

2 

0 

-a 

40/81 

3 

a 

-a 

25/81 

4 

-a 

0 

40/81 

5 

0 

0 

64/81 

6 

a 

0 

40/81 

7 

-a 

a 

25/81 

8 

0 

a 

40/81 

9 

a 

a 

25/81 

a  =  0.7746 


FIGURE  6.5.  SAMPLE  POINTS  AND  WEIGHT  FACTORS  FOR  GAUSSIAN 
QUADRATURE  OF  EQUATION  (6.25) 
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the  terms  from  equation  (6.25)  to  those  given  by  equation  (6,23). 

The  terms  of  the  element  forcing  matrix  are  given  by  the  right 
hand  side  of  equation  (6.13).  Since  the  first  two  terms  in  the  surface 
integral  deal  with  flow  through  the  external  boundaries  of  the  flow 
domain  they  can  be  lumped  together  and  treated  as  one.  These,  in  effect, 
describe  the  boundary  conditions  of  the  problem.  It  will  be  assumed  that 
the  only  impressed  flow  on  the  boundary  of  the  flow  domain  will  be  due 
to  the  presence  of  a  well.  The  domain  can  be  subdivided  in  such  a  way 
that  the  only  element  surfaces  in  contact  with  the  well  are  those  which 
contain  nodes  1,  2,  5,  6,  9,  13,  14,  and  17.  If  the  restriction  is 
placed  on  this  surface  that  it  is  always  rectangular  then  the  integral 
of  the  flow  from  the  well  over  this  surface  can  be  evaluated  once  and 
for  all  to  give  the  distribution  of  this  flow  among  the  nodes  concerned. 
The  terms  of  the  element  forcing  matrix  due  to  well  flow,  then,  are  given 

by 

f j  =  Q  I        f    Nj  |det  J|  di;  dC  (6.26) 

-1   -1 

where  Q  is  the  flow  into  the  element  and  is  assumed  to  be  evenly  dis- 
tributed across  the  surface.  This  integral  is  evaluated  numerically 
using  the  four  sample  points  and  the  weighting  functions  shown  in  Figure 
6.6.  The  results  of  this  integration  are  that  the  inflow  Q  should  be 
distributed  among  the  eight  nodes  as  follows 


fl  =  ^2  =  ^5  =  ^6  =  -Q/"'^  (^-^^a) 

fg  =  -^13  =  fi4  =  fi7  =  Q/3  (6.27b) 
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FIGURE  6.6.  SAMPLE  POINTS  FOR  GAUSSIAN  QUADRATURE  OF  EQUATION  (6.25) 
a  =  0.57735;  WEIGHT  FACTOR  =  1 
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This  somewhat  surprising  result  is  due  to  the  fact  that  a  non-linear 
interpolation  function  is  being  used. 

The  last  term  in  the  integral  on  the  right  hand  side  of  equation 
(6.13)  is  part  of  the  leakage  term  in  the  governing  equation.  It  is 
integrated  over  the  top  surface  of  the  leaky  elements  in  much  the  same 
way  as  for  the  leakage  terms  in  the  element  stiffness  matrix  except 
that  only  four  sampling  points  are  used. 

Computer  Implementation 

The  previous  section  showed  how  the  individual  terms  of  the  element 
matrices  are  computed.  For  each  element  this  results  in  a  20  x  20 
element  stiffness  matrix,  [k],  and  a  1  x  20  forcing  matris,  {f}.  After 
these  element  matrices  are  obtained  they  must  be  inserted  into  the 
correct  locations  in  the  corresponding  global  matrices.  If  n  is  the 
number  of  nodal  points  which  have  been  used  to  subdivide  the  flow  domain, 
then  the  global  matrices  consist  of  an  n  x  n  stiffness  matrix,  [K],  a 
1  X  n  forcing  matrix,  {F},  and  a  1  x  n  matrix  of  nodal  unknowns,  {$}. 
These  matrices,  when  arranged  in  the  form  of  equation  (6.15)  represent 
a  system  of  n  algebraic  equations  in  n  unknowns  which  can  be  solved  by 
any  of  several  methods.  In  choosing  a  method  of  solution  it  is  useful 
to  observe  that  many  of  the  terms  the  stiffness  matrix,  [k],  are  zero. 
In  fact,  if  the  flow  domain  has  been  subdivided  into  elements  in  a 
propitious  manner,  all  of  the  non-zero  terms  will  lie  in  a  relatively 
narrow  band  centered  on  the  diagonal  of  the  matrix.  If  the  method  of 
Gaussian  elimination  is  used  to  solve  these  equations  the  terms  of  the 
stiffness  matrix  which  lie  outside  this  band  need  not  be  dealt  with  at 
all.  This  makes  the  solution  of  the  system  of  equations  by  Gaussian 
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elimination  much  more  efficient  than  solution  by  inversion  of  the  entire 
matrix.  This  advantage  becomes  even  more  pronounced  when  use  is  made  of 
the  symmetric  nature  of  this  matrix.  This  means  that  the  actual  band- 
width to  be  used  consists  only  of  the  diagonal  terms  and  the  terms  on 
one  side  of  the  diagonal  within  the  non-zero  strip.  The  bandwidth  can 
be  minimized  by  subdividing  the  region  and  numbering  the  nodes  in  such 
a  way  that  the  difference  between  the  largest  and  smallest  node  numbers 
in  any  element  is  as  small  as  possible.  For  elements  with  four  nodes  or 
less  computer  programs  are  available  which  automatically  number  the 
nodes  to  give  minimum  bandwidth  (Collins,  1973).  When  the  elements  con- 
tain 20  nodes,  however,  bandwidth  minimization  is  left  to  the  ingenuity 
of  the  programmer. 

Program  Verification 
In  order  to  demonstrate  the  validity  of  the  finite  element  model, 
a  groundwater  flow  problem  for  which  an  analytical  solution  is  available 
was  chosen  for  comparison.  This  is  the  problem  of  axisymmetric  flow 
from  an  injection  well  in  a  bounded  leaky  aquifer,  which  was  described 
in  Chapter  5.  The  solution  is  given  in  terms  of  modified  Bessel  func- 
tions by  equation  (5.38).  If  the  position  of  the  aquifer  boundary  is 
r^  and  the  push  up,  s,  at  that  point  is  zero  then  equation  (5.38)  becomes 


Q, 


f  K  (r   /B) 


(6.28) 


which  is  the  solution  given  by  Jacob  (1946). 

The  analysis  which  led  to  equation  (6.28)  differs  from  the  finite 
element  model  in  that  it  is  one  dimensional  in  radial  coordinates  while 
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the  model  is  three  dimensional  in  Cartesian  coordinates.  However,  when 
the  value  of  the  leakage  coefficient,  B,  is  large  the  flow  is,  in  fact, 
nearly  radial  and  the  two  solutions  should  coincide.  As  B  gets  smaller, 
the  one  dimensional  solution  becomes  less  valid  and  this  comparison  will 
show  the  limits  of  its  accuracy. 

The  shape  of  the  flow  domain  is  shown  in  Figure  6.7.  It  is  divided 
into  seven  elements  with  92  nodes.  Since  the  flow  is  axisymmetric,  the 
sides  and  bottom  of  the  wedge  are  boundaries  across  which  there  is  no 
flow.  Leakage  is  permitted  across  the  top  surface  of  the  wedge  and  there 
is  a  specified  flow  rate  into  the  narrow  end  from  the  well.  There  is  an 
unspecified  amount  of  flow  through  the  outer  end  of  the  wedge  where  the 

drawdown  is  maintained  at  s  =  0.  The  injection  rate  into  the  well  is 

3 
0.1  m  /s  and  since  the  flow  is  axisymmetric  the  wedge  receives  1/12  of 

3 
this,  or  0.00833  m  /s.  Since  the  thickness  of  the  aquifer  is  b  =  10  m 

and  the  hydraulic  conductivity  is  K  =  0.001  m/s  the  one  dimensional  and 

three  dimensional  leakage  coefficients  are  related  by 

B  =  AF  6  =  0.13  (6.29) 

The  solutions  for  push  up,  s,  given  by  the  two  methods  for  several 
values  of  the  leakage  coefficient,  B,  are  shown  in  Figure  6.8.  For  B 
values  of  60  m  and  greater  the  two  solutions  are  indistinguishable  at 
this  scale.  For  the  smaller  values  of  B,  the  leakage  rate  is  higher  and 
the  one  dimensional  solution  gives  lower  values  than  the  computer  model. 
In  these  cases  of  high  leakage  rate,  the  vertical  flow  components  are 
large  enough  to  make  the  pressure  distribution  deviate  noticeably  from 
the  hydrostatic.  The  solutions  plotted  for  these  cases  are  based  on  the 
piezometric  head  at  mid  depth  in  the  aquifer. 
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Freshwater  Lens  in  a  Leaky  Saline  Aquifer 
The  logical  extension  of  the  injection  well  problem  in  a  bounded 
aquifer  is  to  check  the  solution  given  in  Chapter  5  for  an  axi symmetric 
freshwater  lens  in  a  leaky  aquifer.  This  can  be  done  with  the  finite 
element  model  by  successive  repositioning  of  the  interface  and  resolving 
of  the  boundary  value  problem  until  the  interface  boundary  condition  is 
satisfied. 

As  before,  a  30°  segment  of  the  flow  domain  was  modeled.  The  shape 
of  this  segment  is  shown  in  Figure  6.9.  Repositioning  of  the  interface 
is  accomplished  in  two  ways.  First,  the  nodes  are  moved  radially  to 
increase  or  decrease  the  size  of  the  lens.  This  is  done  by  comparing 
pressure  computed  for  the  node  at  the  toe  of  the  interface  with  the 
pressure  that  would  be  required  to  satisfy  the  interface  boundary  condi- 
tion. If  the  computed  pressure  is  too  high,  then  the  distance  to  the 
toe,  r^,  is  increased.  Since  the  radii  of  all  of  the  other  nodes  are 
computed  in  proportion  to  r^  this  increases  the  size  of  the  entire  lens. 
If  the  computed  pressure  is  too  low,  r.  is  reduced. 

The  second  adjustment  of  the  interface  position  is  done  by  changing 
the  elevation  of  all  of  the  interfacial  nodes  except  for  the  nodes  on 
the  interface  toe.  The  revised  elevation  of  each  of  these  nodes  is 
obtained  from  equation  (3.17)  based  on  the  head  computed  for  the  node 
during  the  previous  trial. 

It  was  found  that  the  application  of  both  of  these  adjustments  for 
each  trial  resulted  in  instability  unless  the  initial  guess  of  interface 
position  was  very  good.  However,  when  the  two  methods  of  adjustment  were 
alternated,  the  scheme  converged  to  the  correct  position  in  a  fairly 
regular  manner.  This  is  not  to  say  that  convergence  was  always  rapid. 
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Often  15  to  20  iterations  were  required  before  a  stable  interface 
position  could  be  reached.  Actually,  a  \/ery   good  initial  guess  of 
interface  position  could  have  been  made  using  the  method  of  Chapter  5. 
However,  since  the  purpose  of  this  analysis  was  to  verify  the  accuracy 
of  the  two  methods  it  was  thought  more  advisable  to  arrive  at  the  answer 
by  totally  independent  means. 

Figures  6.10  and  6.11  show  the  comparison  of  interface  locations 
arrived  at  by  the  two  methods  under  two  different  cases  of  injection 
rate  and  leakage  coefficient. 
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CHAPTER  7 
NUMERICAL  SIMULATION  OF  LENS  GROWTH 

General 


Chapters  5  and  6  of  this  report  have  been  concerned  with  locating 
the  boundaries  of  freshwater  lenses  under  steady  flow  conditions.  The 
boundary  between  freshwater  and  saltwater  was  easily  definable  in  such 
cases  because  there  was  no  flow  component  normal  to  it  and,  consequently, 
no  longitudinal  dispersion.  The  plan  of  attack,  then,  was  to  treat  the 
freshwater  lens  as  a  distinct  flow  domain  separated  from  the  resident 
saltwater  by  an  interface.  In  the  present  chapter  the  problem  of 
unsteady  flow  will  be  considered.  Here  there  is  a  large  component  of 
flow  normal  to  the  outer  boundary  of  the  lens  and  longitudinal  dispersion 
becomes  significant.  Since  no  interface  can  be  distinguished  to  separate 
the  injected  fluid  from  the  resident  fluid,  it  is  necessary  to  treat 
this  as  a  single  flow  domain  containing  a  non-homogeneous  fluid. 

The  two  equations  which  describe  the  movement  of  the  freshwater 
lens  are  the  governing  equation  for  groundwater  flow,  given  in  Chapter  3 
as  equation  (3.8),  and  the  convective  dispersion  equation,  given  as 
equation  (3.13).  If  the  injected  freshwater  behaved  as  a  tracer-like 
substance,  the  application  of  these  equations  would  be  relatively 
straightforward.   Problems  involving  the  dispersion  of  tracers  have 
been  solved  both  analytically  and  numerically  by  Hoopes  and  Harleman 
(1967),  Ogata  and  Banks  (1961)  and  others.  The  injected  freshwater, 
however,  differs  from  the  resident  saltwater  in  density  and,  to  a  lesser 
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degree,  in  viscosity.  Consequently,  the  changes  in  fluid  composition 
predicted  by  the  dispersion  equation  have  a  direct  influence  on  the 
parameters  in  the  flow  equation.  In  order  to  make  these  influences 
explicit,  it  is  convenient  to  write  the  flow  equation  in  terms  of 
density,  viscosity  and  pressure  as 

V  .  £|i(VP-  pgvz)  -^  (7.1) 

This  differs  from  equation  (3.8)  in  that  hydraulic  conductivity  is 
expressed  in  terms  of  intrinsic  permeability  and  the  fluid  properties, 
piezometric  head  is  expressed  in  terms  of  pressure  and  elevation,  and 
the  change  in  the  mass  of  stored  fluid  is  expressed  directly. 

The  convective  dispersion  equation  is  expanded  in  a  similar 
fashion  as 

-^   (npc)  -  V-  [pc  ^  (VP  -  pgVz)]  =  V-  (pD)  •  Vc      (7.2) 

These  two  equations  are  coupled  through  the  dependence  of  density  and 
viscosity  on  the  saltwater  concentration  as  illustrated  in  the  section 
on  fluid  properties  in  Chapter  3.  This  set  of  coupled  equations, 
together  with  the  constituative  relations  for  fluid  properties,  must  be 
solved  numerically.  It  happens  that  a  computer  program  is  already 
available  (INTERCOM?,  1976)  which  can  applied  to  this  problem.  The 
three  dimensional  solutions  for  transient  lens  flow  obtained  with  this 
numerical  model  are  valuable  in  themselves  as  well  as  in  providing  a 
check  on  existing  two  dimensional  solutions. 

Description  of  the  INTERCOMP  Model 
The  INTERCOMP  model  was  developed  for  the  U.S.  Geological  Survey 
by  INTERCOMP  Resource  Development  and  Engineering,  Inc.,  of  Houston, 
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Texas.  It  is  a  three-dimensional  finite  difference  model  which  was 
designed  to  predict  the  effects  of  liquid  waste  disposal  by  deep  well 
injection  into  saline  aquifers.  It  is  a  model  of  very   general  appli- 
cability to  flow  in  confined  aquifers  under  a  wide  variety  of  boundary 
conditions.  Basically,  it  consists  of  a  numerical  solution  to  equations 
(7.1)  and  (7.2)  together  with  the  equation  of  energy  transport  by 
convection  and  conduction.  Since  the  influence  of  temperature  is  being 
neglected  in  this  report,  only  that  part  of  the  computer  program 
which  deals  with  equations  (7.1)  and  (7.2)  will  be  used. 

The  difficulties  involved  in  applying  finite  difference  approxima- 
tions to  equation  (7.2)  have  been  well  documented  (von  Neumann  and 
Richtmyer,  1950;  Stone  and  Brian,  1963).  Basically,  the  problem  is 
that  if  the  finite  difference  approximation  to  either  the  time  deriva- 
tive or  the  convective  spatial  derivative  is  correct  only  to  the  first 
order,  then  the  resulting  second  order  truncation  error  has  the  same 
form  as  the  dispersive  term  of  the  equation.  This  error  is  called 
numerical  dispersion  and  its  magnitude  can  be  such  that  it  masks  the 
actual  physical  dispersion  being  modeled.  If  a  second  order  correct 
differencing  scheme  is  used,  the  truncation  error  no  longer  has  the  same 
form  as  the  dispersive  term  but  a  restriction  must  be  imposed  on  the 
time  and  distance  steps  used  in  order  to  ensure  the  stability  of  the 
solution.  The  actual  form  of  this  restriction  depends  on  the  differ- 
encing scheme  being  used. 

The  INTERCOMP  model  provides  users  with  a  choice  of  differencing 
schemes.  For  the  computational  runs  mentioned  in  this  report,  the  time 
derivatives  were  represented  by  the  Crank-Nicholson  approximation  and 


130 


the  spatial  derivatives  by  the  central  difference  approximation.  Since 
both  of  these  differencing  methods  are  correct  to  the  second  order,  their 
application  results  in  no  numerical  dispersion.  However,  to  ensure 
stability,  the  time  and  distance  steps  should  be  chosen  to  satisfy  the 
following  requirement: 


^    -1-0  (7.3) 


Good  results  can  often  be  obtained  even  if  this  criterion  is  exceeded, 
but  they  are  not  ensured. 

The  Crank-Nicholson  differencing  method  is  an  implicit  one  because 
it  makes  use  of  unknown  values  of  the  dependent  variable.  This  results 
in  large  system  of  simultaneous  algebraic  equations  which  must  be  solved 
at  each  time  step.  The  model  can  solve  these  equations  either  by 
reduced  bandwidth  Gaussian  elimination  or  by  successive  overrelaxation. 
In  this  report,  the  latter  method  was  used. 

The  INTERCOMP  model  has  the  capability  to  solve  either  three 
dimensional  problems  in  cartesian  coordinates  or  axisymmetric  problems 
in  cylindrical  coordinates.  Since  the  problems  to  be  solved  here  deal 
with  flow  from  a  single  well,  the  cylindrical  coordinate  system  is  used. 
For  axisymmetric  flow  in  a  confined  aquifer  the  boundary  conditions 
consist  of  no-flow  boundaries  at  the  top  and  bottom  of  the  aquifer,  a 
prescribed  rate  of  inflow  at  the  well,  and  a  constant  head  boundary  at 
infinity.  The  boundary  condition  at  infinity  is  simulated  by  calcula- 
ting the  rate  of  efflux  from  the  peripheral  grid  blocks  based  on  the 
rate  of  change  of  pressure,  the  compressibility  of  the  fluid  and  porous 
medium,  and  the  total  pore  volume  of  the  modeled  portion  of  the  aquifer. 
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More  sophisticated  methods  for  simulating  infinite  aquifers  are  available 
(Carter  and  Tracy,  1960),  but  since  pressure  distribution  is  of  more 
interest  than  absolute  pressure  in  the  analysis  of  freshwater  lenses, 
this  was  considered  sufficient. 

Input  data  for  this  model  includes  specification  of  the  grid  size, 
time  step,  and  method  of  solution  as  well  as  the  properties  of  the  porous 
medium,  resident  fluid,  and  injected  fluid.  The  finite  difference  grid 
used  in  the  various  computer  runs  is  shown  in  Figure  7.1.  The  grid 
points  shown  in  this  figure  are  computational  nodes  which  represent  the 
center  of  the  grid  blocks  into  which  the  flow  region  is  divided.  The 
pressure  and  concentration  at  each  of  these  points  is  computed  at  ewery 
time  step.  The  modeled  region  is  subdivided  into  10  horizontal  strips 
of  equal  Az  and  25  vertical  strips  grouped  in  three  zones  of  equal 
Mn{r).     The  time  step  used  with  this  grid  was  t  =  0.1  day. 

A  series  of  simulations  using  different  injection  rates  and  aquifer 
configurations  was  run  in  order  to  show  the  effects  of  the  different 
flow  parameters  on  the  shape  of  the  freshwater  lens.  Those  properties 
of  the  aquifer  and  the  resident  and  injected  fluids  which  were  not 
varied  are  listed  in  Table  7.1. 

Application  of  the  numerical  model  results  in  a  list  of  the  pres- 
sures and  concentrations  for  all  of  the  grid  points  at  selected  time 
steps.  In  these  results  pure  freshwater  is  assigned  a  concentration  of 
one,  and  undiluted  saltwater  a  concentration  of  zero.  From  the  model's 
viewpoint,  the  only  difference  between  these  two  fluids  is  their  density. 

Model  Verification 
The  accuracy  of  the  INTERCOMP  model,  when  used  for  its  intended 
purpose,  was  verified  for  several  test  cases  by  its  developers 
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(INTERCOMP,  1976).  It  was  also  tested  by  Papadopulos  and  Larson  (1978) 
against  field  data  for  a  heated  water  storage  well  (Molz  et  al. )  which 
it  matched  satisfactorily.  In  order  to  verify  the  model  with  the  grid 
shown  in  Figure  7.1  and  the  input  data  of  Table  7.1,  a  test  case  was 
used  in  which  the  density  of  the  injected  water  was  the  same  as  that  of 
the  resident  saltwater.  Thus  the  problem  was  reduced  to  the  case  of 
constant  injection  of  a  tracer-like  substance  for  which  approximate 
analytical  solutions  are  available.  The  approximate  solution  given  by 
Hoopes  and  Harleman  (1965)  for  injection  of  a  tracer  with  concentration 
c  =  1  into  an  aquifer  with  initial  concentration  c  =  0  everywhere  is 


c  -  1/2  erfc 


r72  -  t' 
./4/3  r'3_ 


(7.4) 


where 


r'  =  r/a 


2irbna 


H 


erfc 


t  = 


longitudinal  dispersivity 

the  complimentary  error  function 

time  since  the  start  of  injection 


A  comparison  of  the  concentration  profile  predicted  by  equation 
(7.4)  with  that  predicted  by  the  INTERCOMP  model  for  the  same  flow  condi- 
tions is  shown  in  Figure  7.2.  It  should  be  mentioned  that  equation  (7.4) 
is  only  an  approximate  solution  because  it  was  obtained  by  setting  the 
initial  rate  of  change  of  concentration,  9c/9t,  equal  to  zero  at  the 
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well,  which  was  represented  as  a  line  source.  Because  of  this,  it 
predicts  the  presence  of  a  finite  amount  of  tracer  in  the  aquifer  at 
time  t  =  0.  This  initial  error  remains  constant  in  the  approximate 
solution  as  time  progresses  and  results  in  the  predicted  edge  of  the 
tracer  cloud  being  farther  from  the  well  than  it  should  be.  The  exact 
solution  to  the  problem  was  given  by  Ogata  (1958)  but,  unfortunately, 
it  contains  an  integral  which  cannot  be  explicitly  evaluated.  Hoopes 
and  Harleman  evaluated  Ogata's  solution  numerically  for  a  few  specific 
cases  and  compared  it  to  their  approximate  solution.  The  results  they 
obtained  in  this  way  look  identical  to  the  comparison  shown  in  Figure 
7.2. 

Lens  Formation  in  Aquifers  With  Weak  Vertical  Flow 
The  experimental  results  presented  in  Chapter  4  indicated  that  the 
thickness  of  the  saline  aquifer  may  have  a  significant  influence  on  the 
shape  of  the  freshwater  lens.  Consequently,  the  initial  simulations 
were  done  with  a  relatively  thin  aquifer  in  which  vertical  flow  compo- 
nents would  have  little  chance  to  develop.  It  was  expected  that  under 
these  conditions  the  formation  of  the  freshwater  lens  would  proceed  as 
predicted  by  previous  investigators  who  neglected  vertical  flow  compo- 
nents entirely  (Kimbler,  Kazman  and  Whitehead,  1975).  Figure  7.3  shows 
the  results  of  numerical  simulation  of  a  freshwater  lens  in  a  10-meter 
thick  aquifer  after  40  days  of  injection.  The  left  border  of  the 
figure  can  be  taken  as  the  location  of  the  well.  The  mixing  zone  is 
represented  as  that  part  of  the  flow  region  between  the  90%  and  10% 
concentration  lines,  although  the  selection  of  these  concentrations 
was  arbitrary.  The  50%  concentration  line  can  be  considered  to  represent 
the  position  of  an  imaginary  interface  which  would  exist  if  there  were 
no  dispersion.  Note  that  the  lines  of  equal  concentration  are  straight 
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and  that  after  40  days  of  injection  their  departure  from  the  vertical 

is  still  very  slight.  In  this  case,  the  amount  of  freshwater  which  could 

be  recovered  from  the  lens  is  limited  primarily  by  dispersion. 

The  thinness  of  the  saline  aquifer  suppresses  gravity  segregation 
of  the  fluids  in  two  ways.  First,  the  pressure  differential  caused  by 
the  difference  in  fluid  density  is  very   small  in  a  thin  aquifer  compared 
to  the  pressure  gradients  created  by  the  injection  well.  Second,  the 
opportunity  for  vertical  flow  components  to  develop  is  limited  by  the 
small  vertical  dimensions  of  the  aquifer.  If  the  aquifer  were  to  be 
made  thicker  and  the  vertical  hydraulic  conductivity  were  reduced,  then 
only  the  second  of  these  suppressing  influences  would  be  in  effect. 
This  was  done  in  the  numerical  simulation  illustrated  by  Figure  7.4. 
The  aquifer  in  this  case  was  50  meters  thick  but  the  vertical  hydraulic 
conductivity  was  reduced  by  an  order  of  magnitude.  The  injection  rate 
was  also  increased  so  that  the  horizontal  velocities  would  be  similar 
to  those  in  the  previous  run.  In  this  case,  the  lines  of  equal  concen- 
tration are  slightly  curved  but  their  rate  of  rotation  toward  the 
horizontal  is  not  significantly  greater. 

Lens  Formation  in  Aquifers  With  Strong  Vertical  Flow 
It  would  be  interesting,  now,  to  see  the  effect  of  the  vertical 
flow  components  on  the  formation  of  a  lens  in  a  thick  aquifer  in  which 
they  are  not  suppressed.  Figure  7.5  shows  the  results  of  a  simulation 
which  was  identical  to  the  previous  one  in  every   way  except  that  the 
vertical  hydraulic  conductivity  was  made  equal  to  the  horizontal.  The 
difference  between  this  figure  and  Figure  7.4  indicates  that  the  vertical 
flow  components  in  a  thick  aquifer  can  be  quite  strong.  Their  effect 
is  to  increase  the  rate  of  gravitational  segregation  and  to  make  the 
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mixing  zone  much  wider  at  the  top  than  at  the  bottom.  The  reason  for 
this  latter  phenomenon  is  revealed  in  Figure  7.6  which  illustrates  the 
movement  of  the  50%  concentration  line  during  the  40  days  of  injection. 
The  lower  part  of  this  lens  remained  almost  stationary  after  the  10th 
day.  Evidently,  the  flow  in  this  region  was  primarily  tangential  to  the 
equal  concentration  lines  so  that  only  lateral  dispersion  was  active. 
The  longitudinal  dispersion  associated  with  the  growth  of  the  lens  was 
limited  to  its  upper  part.  From  the  viewpoint  of  freshwater  storage, 
the  conditions  illustrated  in  Figure  7.6  are  even  worse  than  they  appear. 
Because  of  the  cylindrical  geometry  of  the  lens,  the  volume  of  water 
stored  is  proportional  to  the  square  of  the  lens  radius.  The  greater 
part  of  the  injected  freshwater  is  stored  beyond  the  toe  of  the  lens  and 
cannot  be  recovered  undiluted  by  pumping  from  the  well. 

Dimensionless  Representation  of  Lens  Growth 
It  is  evident  from  the  numerical  simulations  that  storage  of  fresh- 
water in  lenses  is  much  more  effective  in  thin  saline  aquifers  than  in 
thick  ones.  Since  the  concept  of  thickness  is  relative,  it  is  convenient 
to  present  the  numerical  results  in  a  dimensionless  form  which  provides 
a  basis  for  comparing  one  aquifer  to  another.  This  is  done  by  defining 
a  dimensionless  time,  t',  an  average  radius  r,  and  an  aquifer  coeffi- 
cient, X,  as 

t'  =  Qt/b^  (7.5a) 

A  =  bh/Q  (7.5c) 

The  results  of  a  series  of  numerical  simulations  are  presented  in  terms 
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of  these  variables  in  Figure  7.7.  The  absissa  of  this  chart  gives  values 
of  v^/v,   where  r^  is  the  radius  of  the  toe  of  the  50%  concentration  line 
and  r,  as  given  by  equation  (7.5b),  would  be  the  radius  of  this  line  if 
there  were  no  gravitational  segregation. 

Figure  7.7  indicates  that  aquifers  with  small  values  of  X  are  more 
favorable  to  the  storage  of  freshwater.  The  applicability  of  this 
figure  for  predictive  purposes,  however,  is  limited  to  homogeneous  and 
isotropic  aquifers  with  a  longitudinal  dispersivity  of  6.1  meters  and  a 
Ghyben-Herzberg  ratio  of  6  =  40. 

The  Influence  of  Dispersion  on  Lens  Shape 
Previous  investigators  have  found  that  the  rate  of  gravitational 
segregation  of  fluids  of  different  density  is  retarded  by  dispersion 
across  the  interface  between  them  (Gardner,  Downie  and  Wyllie,  1962; 
Esmail  and  Kimbler,  1967).  This  effect  was  included  in  Whitehead's 
analysis  of  freshwater  lenses  in  thin  aquifers  (Whitehead,  1974). 
Figure  7.8  shows  a  comparison  of  two  freshwater  lenses  after  30  days  of 
injection  under  conditions  that  were  identical  except  for  the  disper- 
sivity of  the  porous  medium.  The  difference  in  the  amount  of  dispersion 
has  an  obvious  effect  on  the  shape  of  the  lens.  The  rate  of  rotation  of 
the  equal  concentration  lines  toward  the  horizontal  is  reduced  by 
doubling  the  dispersivity.  For  very  extensive  and  thin  lenses  and  long 
term  storage  this  could  have  a  significant  effect  on  the  recovery 
efficiency.  For  the  case  shown,  however,  the  reduction  in  gravity 
segregation  was  balanced  by  an  increase  in  the  width  of  the  mixing  zone 
so  that  the  amount  of  freshwater  which  could  be  recovered  was  not  greatly 
changed. 
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FIGURE  7.7.  A  DIMENSIONLESS  REPRESENTATION  OF  GRAVITATIONAL 
■   SEGREGATION  DURING  LENS  GROWTH 
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Changes  in  Lens  Shape  During  Recovery 
The  production  of  freshwater  from  a  lens  is  accompanied  by  the 
same  irreversible  processes  of  dispersion  and  gravity  segregation  that 
were  present  during  its  injection.  The  only  difference  is  that  during 
production  the  lens  boundaries  are  flowing  in  the  direction  of  increas- 
ing pressure  gradient  so  that  the  rate  of  movement  of  the  boundaries 
increases  rapidly  as  the  well  is  approached.  Figures  7.9  and  7.10 
illustrate  the  collapse  of  freshwater  lenses  during  production.  The 
aquifer  of  Figure  7.9  has  a  X   value  of  2.5  so  conditions  are  relatively 

favorable  for  the  storage  of  freshwater.  The  lens  in  this  figure  was 

3 
formed  by  injection  of  freshwater  at  a  rate  of  0.1  m  /s  for  30  days. 

Removal  of  freshwater  from  the  lens  was  begun  on  the  31st  day  at  the 

same  pumping  rate.  By  the  54th  day,  the  50%  concentration  line  had 

broken  through  to  the  well.  The  amount  of  water  produced  at  that  time 

was  74%  of  the  amount  that  had  been  injected  although  the  produced 

water  was  not  100%  fresh. 

Figure  7.10  illustrates  the  same  sequence  of  events  in  an  aquifer 

having  a  X   value  of  22.05.  In  this  case  the  50%  line  reached  the  well 

by  the  38th  day  and  the  recovery  efficiency  at  this  time  was  27%. 
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CHAPTER  8 
CONCLUSIONS 

1.  The  equation  for  predicting  the  rate  of  rotation  of  an 
initially  vertical  interface  between  two  fluids  of  different  density 
given  by  Gardner,  Downie  and  Kendall  (1962)  matches  the  experimental 
data  from  the  narrow  Hele-Shaw  model  quite  well.  The  shape  of  this 
interface  as  it  rotates  toward  the  horizontal  is  nearly  straight  as 
assumed  in  the  derivation  of  the  equation.  When  the  depth  of  the  flow 
region  is  doubled,  however,  the  moving  interface  becomes  more  curved 
and  the  rotation  rate  is  greater  than  predicted.  This  is  probably 
caused  by  the  increased  opportunity  for  vertical  flow  in  the  thicker 
model  aquifer. 

2.  Saltwater  intrusion  in  coastal  aquifers  can  be  alleviated  on 
a  local  scale  by  the  creation  of  a  freshwater  lens  or  a  series  of  such 
lenses.  Such  measures  could  be  especially  useful  in  highly  developed 
coastal  zones  where  the  natural  recharge  of  the  coastal  aquifer  is 
reduced  by  land  drainage  practices  and  the  construction  of  impermeable 
land  surfaces.  An  injection  well  scheme  of  this  type  could  be 
evaluated  and  designed  using  the  Girinskii  potential  method  presented 
in  Chapter  5. 

3.  It  seldom  happens  in  nature  that  an  aquifer  is  confined  by 
completely  impervious  layers.  Many  aquifers  that  might  be  considered 
as  possible  sites  for  freshwater  lens  storage  projects  are  actually 
more  or  less  leaky.  The  amount  of  leakage  that  would  be  caused  by 
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the  injection  of  freshwater  has  an  important  influence  on  the  size 
of  freshwater  lens  that  can  be  formed.  The  limiting  size  of  lens  that 
can  be  formed  under  a  given  set  of  flow  conditions  in  a  leaky  aquifer 
can  be  determined  using  the  design  chart  developed  in  Chapter  5. 

4.  The  classical  one-dimensional  solution  for  flow  in  leaky 
aquifers  which  was  developed  by  Jacob  (1946)  gives  essentially  the  same 
results  as  the  three-dimensional  numerical  solution  as  long  as  the 
ratio  of  B/b  is  greater  than  10. 

5.  The  three-dimensional  finite  element  model  developed  in  Chapter 
6  gives  yery   stable  and  reliable  solutions  to  the  groundwater  flow 
equation  within  the  curved  boundaries  typical  of  a  freshwater  lens. 

The  extra  programming  effort  required  by  the  use  of  quadratic  inter- 
polation functions  and  the  20  node  elements  is  rewarded  by  the  larger 
element  size  and  more  flexible  element  shape  it  makes  possible. 

6.  The  location  of  a  stationary  interface  can  be  found  by 
iterative  solution  of  the  boundary  value  problem  within  the  freshwater 
lens  using  the  three-dimensional  finite  element  model.  If  the  resident 
saltwater  is  also  in  motion,  the  same  technique  can  be  used  but  two 
boundary  value  problems  must  be  solved  at  each  iteration.  Unless  the 
initial  guess  of  the  location  of  the  interface  is  very   good,  a  large 
number  of  iterations  may  be  needed. 

7.  The  INTERCOMP  finite  difference  model  is  very  useful  in 
tracing  the  movement  of  a  transient  freshwater  lens.  This  model  is 
capable  of  handling  non-homogeneous  and  non-isotropic  aquifers.  If 
the  temperature  variations  within  the  aquifer  are  substantial,  it 
would  be  desirable  to  include  their  effects  also.  This  can  be  done 
with  the  INTERCOMP  model.  It  would  certainly  be  worthwhile  for  anyone 
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planning  a  large  scale  lens  storage  project  to  check  its  feasibility 
with  this  model  first. 

8.  The  vertical  components  of  flow  within  a  freshwater  lens  can 
have  a  ^ery   powerful  influence  on  the  shape  of  the  lens.  Previous 
analyses  have  assumed  that  vertical  flow  components  are  negligible. 
This  must  be  true  if  the  lens  is  to  be  useful  for  freshwater  storage, 
but,  in  general,  it  is  not  true.  In  thick  aquifers  with  good  vertical 
hydraulic  conductivity,  the  vertical  components  of  flow  have  a 
predominant  influence  and  under  these  conditions  the  feasibility  of 
storage  in  freshwater  lenses  is  highly  questionable. 

9.  The  concept  of  storing  water  in  freshwater  lenses  for  future 
use  is  based  on  an  idealization  of  groundwater  flow.  When  the  natural 
conditions  conform  closely  to  these  idealized  conditions,  reasonably 
high  recovery  efficiency  can  be  achieved.  With  each  new  departure 
from  the  ideal  conditions,  however,  the  recovery  efficiency  is  reduced 
and  the  scheme  rapidly  becomes  unfeasible. 


■APPENDIX 
FINITE  ELEMENT  PROGRAM  DOCUMENTATION 

Input  Data 
The  required  input  data  for  this  program  consists  of  an  initial 
data  list  and  a  repeated  data  list.  The  initial  data  list  is  read  by 
the  main  program  and  consists  of  the  following  variables: 

IBWTH  =  Bandwith  of  the  discretization  scheme 

2 
B2    =  3  ,  the  modified  leakage  coefficient 

HO    =  Piezometric  head  in  the  aquifer  above  the  leaky  layer 

NUNODS  =  Number  of  nodes 

NUMELS  =  Number  of  elements 

YKY    =  K 
lU         -   K^^ 

RT    =  Radius  of  the  interface  toe 
DEPTH  =  Thickness  of  the  saline  aquifer 
QQ    =  Rate  of  injection 

M     =  Number  of  the  first  element  with  an  adjustable  boundary 
The  repeated  data  list  is  read  in  subroutine  LDCAT  and  consists 
of  the  following  variables: 

R     =  Ratio  of  node  radius  to  RT 
TH(1)  =  Elevation  of  the  top  tier  of  nodes 
TH(2)  =  Elevation  of  the  middle  tier  of  nodes 
TH(3)  =  Elevation  of  the  bottom  tier  of  nodes 
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Flow  Chart 
The  finite  element  model  consists  of  a  main  program  and  seven 
subroutines.  The  primary  function  of  the  main  program  is  to  control 
the  subroutines  in  which  most  of  the  actual  computations  are 
accomplished. 


(^  start) 


READ  INITIAL 
DATA 


ITS  =  0 


ITT  =  0 


CALL  SUBROUTINE  LOCAT 


INITIALIZE  GLOBAL 
MATRICES  TO  ZERO 


WRITE  INITIAL  DATA,  NODAL 
COORDINATES  AND  BOUNDARY  VALUES 


ITS  =  ITS  +  1 


IE  =  1 
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INITIALIZE  ELEMENT  MATRICES 
TO  ZERO 


FIND  COORDINATES  OF  ALL 
ELEMENT  NODES 


CALL  SUBROUTINE  ISOVER 


■YES 


IE  =  IE  +  1 


NO 


CALL  SUBROUTINE  LEAKY 


ASSEMBLE  ELEMENT  MATRICES  INTO  THE 
GLOBAL  MATRICES 


■NO- 


YES 


APPLY  BOUNDARY  CONDITIONS 
AT  EVERY  NODE 


CALL  SUBROUTINE  BAND 


CALCULATE  PRESSURE  AT 
EVERY  NODE 
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0 


WRITE  NODAL  COORDINATES,  HEADS 
AND  PRESSURES 


YES 


ITS  =  ITS  +  1 


YES 


CALL  SUBROUTINE  AJUSE 


ITT  =  ITT  -  1 


NO 


CALL  SUBROUTINE  AJUST 


ITT  =  ITT  +  1 


© 
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(subroutine  ISOVER^ 


ISPOT  =  1 


SELECT  THE  NEXT  SAMPLE 
POINT  AND  WEIGHTING  FACTOR 


COMPUTE  THE  PARTIAL  DERIVATIVES  OF  THE 
SHAPE  FUNCTIONS  AT  THE  SAMPLE  POINT 


COMPUTE  THE  TERMS  OF  THE 
JACOBIAN  MATRIX 


_L 


COMPUTE  THE  ABSOLUTE  VALUE  OF 
THE  JACOBIAN  DETERMINANT 


CALL  SUBROUTINE  INVERT 


MULTIPLY  THE  INVERSE  JACOBIAN 
MATRIX  BY  THE  PARTIAL  DERIVATIVES 
OF  THE  SHAPE  FUNCTIONS 


DETERMINE  THE  VALUE  OF  THE 
INTEGRAND  AT  THE  SAMPLE  POINT 
AND  ADD  IT  TO  THE  ELEMENT 
STIFFNESS  MATRIX 


(return) 
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(subroutine  leaky) 


INITIALIZE  TEMPORARY  ELEMENT 
MATRICES  TO  ZERO 


I SPOT  =  1 


CHOOSE  NEXT  SAMPLE  POINT 
AND  WEIGHING  FACTOR 


CALCULATE  THE  VALUE  OF  THE 
SHAPE  FUNCTIONS  AT  THE  SAMPLE  POINTS 


COMPUTE  THE  ABSOLUTE  VALUE  OF  THE 
JACOBIAN  DETERMINANT  AT  THE 
SAMPLE  POINT 


DETERMINE  THE  VALUES  OF  THE 
INTEGRANDS  AT  THIS  POINT  AND 
ADD  THEM  TO  THE  ELEMENT  MATRICES 


ISPOT  =  ISPOT  +  1 


■NO- 


c 


YES 


RETURN 


) 
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(subroutine  ajust) 


COMPUTE  THE  PRESSURE 

ON  THE  FRESHWATER  SIDE 

OF  THE  INTERFACE  TOE,  PF 


COMPUTE  THE  PRESSURE  ON 
THE  SALTWATER  SIDE  OF 
THE  INTERFACE  TOE,  PS 


COMPUTE  THE  AVERAGE  HORIZONTAL 
PRESSURE  GRADIENT  IN  THE 
OUTERMOST  ELEMENT  OF  THE 

FRESHWATER  LENS 


RELOCATE  THE  TOE  OF  THE  INTERFACE 
BY  LINEAR  EXTRAPOLATION  TO 
MAKE  PS  =  PF 


RELOCATE  ALL  OF  THE  REMAINING 
NODES  IN  PROPORTION  TO  THE  SHIFT 
OF  THE  INTERFACE  TOE 


(return) 
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(subroutine  ajuse) 


IE  =  M 


I 


FIND  NODES  3  AND  4  OF 

ELEMENT  IE  AND  SHIFT 

THEM  VERTICALLY  TO  SATISFY 

THE  INTERFACE  BOUNDARY 

CONDITIONS 


± 


FIND  NODES  10,  11,  AND  12 
OF  ELEMENT  IE  AND  SHIFT  THEM 

VERTICALLY  TO  SATISFY  THE 
INTERFACE  BOUNDARY  CONDITIONS 


FIND  NODES  15  AND  16  OF 
ELEMENT  IE  AND  RELOCATE  THEM 
HALF  WAY  BETWEEN  THE  NODES 
ABOVE  AND  BELOW  THEM 
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(subroutine  locat) 


THIS  SUBROUTINE  TRANSFORMS  THE  REPEATED 

INPUT  DATA  FROM  THE  CYLINDRICAL  COORDINATE 

SYSTEM  INTO  THE  CARTESIAN  COORDINATE  SYSTEM 


(subroutine  invert) 


THIS  subroutine  INVERTS  THE  JACOBIAN 
MATRIX  BY  THE  STANDARD  GAUSS-JORDAN  METHOD 


(subroutine  band) 


1 

THIS  SUBROUTINE 

SOLVES 

THE  SYSTEM 

OF 

SIMULTANEOUS  ALGEBRAIC 

EQUATIONS 

BY 

REDUCED  BANDWITH 

GAUSSIAN  ELIMINATION 

161 


i: 

1  ■/ 
-.3 


■'  r.  -i  D  '  7 ,  '.  t- 1  "  ■■    "  .C   ■-■  0  •  i'--'.  ,\  ■  ■■  t-  r' ,  ; :;  ' 


'JO    TO    iO 

3"  -0.  X'^y.  VAV,  :kz 


•.P  ITZ'.s,  2C2  = 


r  ! ,  SN">  (  I  ;  .  J'  !  } 


t  ;  -'J     J 


~0    1  =  '. 

'J  :i.  ; 


:alL     rZCVa"  :  x;-..',  -^'/Y-  Z':Z,  '-'CL.'ME) 

>       _£r^  ■■■;?.         i    '         -'J      TH     i5 

:al_    uEAi^v.jr    >-o,  "LA  . 


r 


:'0  =0   1  =  1,  ■■.!.  v-r^ 


■ --Ai 1 . J) 


1=1    0! 


-1     :■     iO    TIT    30 
;  t  1  ; 


■  -  ..I'NGCS; 


I ' ,  F  ■■■  n  .  >= 


"■J ; 


;r.    ~o    ^cQ 


z.  J 


■'  1  ■* 


■--LL    A^',S' '-J'-'-i:.-::.  .■'T-.    !,   .',  CEPTH.  F.  N,  ,~N.  "HN'. 
;T-=!TT.1 

'Z-L^    A.I.  S^r  .  "J,  ■-,  ■j,_:v-£i_3,  ;     rcpTH,-. 


1  ■■  ~-^:3-7-    =,  ;j.  ■:-(,  ^HhLNOOS    =.   [A,   r^.   ■ 
1<     ^f-2c.    =    ui:^    .1.  ^x,  AHrO    =-Fl.:     o,  =:;c.=^ 


162 


35 
96 


-3 

5o 


91 


"a 
-5 


I    :    J 

:  1- 

i  1^ 

I  15 
'  i? 

iZ: 


;o::5    F^rrATi.tX     .-,kO0  =  .€1O    4.  5<.  3HC'  =  .   T-1  1 

2030    ;-Cr5r^AT:.  IHO.  "HcLEr^'^.r '7    5r.,   ■iHT''  =  i=:.  ;x    4h^-.Aa.  ;X.  :t'r-  \  I  J 

IrCHJ  5  -  7         3         ^       i;       11        !2       1~       :4       IT       lo. 

215H    17       1-3       :?       SO/ 
rO?5    P"C''"'iT  '  1 '^.  "C:-' 

IcJiH     .  ,  ... 

^■-■40    -';FM4r' -I'- .  ■-■•  ■"■■  ■    ■  ■- ■  "■■•',  ^'l.  d'<  ,  SOT") 
rCSO    "CRMAKlriC.  -i^MODE.  i'^.  5H       X       .  SJ  .   rh'       V        .  5-'/  T,-'       Z       •  "'. 

i  10"  3^■'-  ,  ;<.  iCH  a 

ZC?5    i^jP'<AT',  l.^.  .tC^  

IZ'jt^ .      .  ,         ) 

2CtO    -QRI'iT;  :,<•  14.  'Sr^lC    -i,  3.-'    F10    7.  5x,f;o    ~) 
SCO    -Oi'fAT .  i.HO,  ;  5H    '0~-i';.    VCL'.'ME    =,EI6    r  ■  .'    'lO''. , 

1 2cH-*>'*-**'**'*-*     ST  ThFi'l^Sc     ^lAfrRI.X     -»■^-^■«'■f»■*-f.•*  *  ; 

I'jSo  Pc~r-A7' I :ei 2  ?■  ■ 

;OSS    rG.-f"^7  ■.  Iril  ,    I'l-i    ""'uTAL    LEAX^'    AREA    =.£!=     5  .•'  ' ,   •'•C.<, 

1  ""7M*  vt*-**-*-*-!!  <      -7Pv^r     .^cT^IX     -♦*-^-tf  •♦■*-*■*-*      ) 
703to    FuK.''A7  ■.  IHO.  -H    ^.CDE    NO.  .  '.OX     ImX,  19.<,  IhY.  •9;<',  IHZ,  It'X    ^HhE.AC  •;.<. 

15l-'=PE.55iJRE..'7v.  ';|.J .  .5,,   ;0K •  10<     lOH . 

210)C.  1C:h .9<.  IIH .=X.  l-M . 

20=0    ~CF:MAr(  4X,  I-l,  Px,  "10.  5,  IGX,  r;0.  5.  :0X,  FIO,  5.  =■<■  S:  .     4    ^X    t.12    bi 
wcr;     3'.jD 


H'JSRGlJ'^  :WE    L£AK  /  '  32,  MO.  '^LA  : 
."  :  .-^ENS  :  G^    E  A  ■:  7'; .  r  ?  :  .  ,^  j  .;  7 ,  ;  ; 
CCM^CW    HA,  Hr  .  S.'<,  EV.  EZ 
30     1     1=1, =0 
£F' I  !-0.  0 
1     i.N'  I  '-0     3 
;  7  =  07=1 

ALF^A  =  0    ^'''i-'iibJ 
C0=75. 0/El     0 

:r<-JO.  0'81.  0 
<  :=-AL?;-^A 

,;=-.iL  =  HA 

'.?  ij     T  C     1  'J 

7  <:=o,  0 

5C    70     10 
3    ,<1=AL?HA 

■;o  to  10 

.•.=C-i 

;.j   -3   10 

7  .X ;  =0  0 

U=s4    C.31    0 

■;o  *Q  -0 

0     a=-ALr^A 


EY' 


2.M;X0) 


-0 
!4{ 


1 43 

144 


.19 
L4=: 


.  31 
52 


v7  TO  :o 
"    '^-.lPma 

■JC    TO    10 

3  xi=o  ■:• 

>n    70    10 

^      ,(  T  =AL,PH.A 

w=0  J 
■•==1     0+rr 


X  It" 


■  •.  ; 


-2=2    oixi ■ 

3;=(  X,7-w  i  ;-*v;.-.!,^     ,-, 

- -5=  (  '.(2-1  :.'•■'  I^  '4     0 

.Tg=  I  x5-y  !  J  •  »  f  ?    -1    0 

F  •.  7- -i  liY'II' 

°  13='  :    0-' I  ■♦'•';  ■'  •  .7.  C 

'•'  '.  '/--  I  I'VIr 

=20=- -15 

3.'-.  1 ,  \  :  -'^ ':■>€.-('  ":    -^?6*Fx  I  3  ;  ■^:■^7♦=;x  ■,  •-  >  -.--V8*E.'  ,  5 
;  j-o  1  7  j»rx  ;  1  7  ."  -w  '_  ;.».7Y  ?  ;  c  ■  4..^  ;  ?.»u  Y  "'^  o  \  j.t;^v;'-i£Y  . 

f  J  ;  1 ,  2  '  =P  5*i  '  ■  5  ■  ■•'^  i  *E  Y  '  -  ,  -^iV  7fl-£  ^  I  "  ^  -^'S^E'' 
I  .:.&  I  7<£  /  ,   ;  7  ;  -;  !  ;  -^i"  ,'  r  ;^5  )  tO  j  .-;»c;v  ,  ;  9  ;,  ->-R20''£''  \  20 

-  7=7    o**-,  ; 

9  ?  =  •  Y  2  -^  <  I "  *  < :  ri ,  4  3 
•"i=(Y2-  > :  )-»'■  i  =  /4  ,:■ 
.=  -=■  vz-^Ki  ;■*'  ;.^.-  4  ■: 
=3=  1  »2-< :  .'■♦X  :>i.  4  J 
=  :■'=(  1   c  -vi*x:  ■>■;-.€ 


IbJ 


I  r  7 
1;3 


lc3 


lc.3 
'-  5-' 

■Zi 


133 
1  ^Jl 


':-!_*-' '  '-  ■■ '  ^^-  '■  5^r  •  ( ;  g  ,  ^?  1  ^^ev  '  I  ?  i  v^2''.*H>  ■  ?o  •■ 

^  ^y  '  5  !  =  -  I  ;  -  »  :•  ;  M  *  .    r  ;  .  y  ;  I,  .      jT  ;    ■  4     ,3 

SN ;  c ;  ==  '■  I  -  ■*  ^  I  :'i  *  •  :<  ■  - ,  1 1 1 '  -, ";  '•  'a  ,-, 
3M(7-.  =  >  ::■»,:=*,  YT. ,■-_;■  Ai/4  rS 

5N.  Birr-,- jv,;,v-o*,  <  l-'I-'-:'".-.  -  a"o 
=N'  (  1  '^  )  =  '  I;-!»  '  1    0-  t  ;  -■'X  ;  '  ,  :>    .'i 
fi- 1  :  3;  =>:?■!•■  1     ,';-vT.>v7  i  '5    ' 

frj  1  ^0  ■  =  <  I  i^s-  ■;  1    :'i —  ;  *Y  i  r  /  ri   n 
LO    10    1  =  1 ,  rr, 
CO    SO    .  =  ■. ,  r.o 

;.:j   :.M    _[  =  i   a 


ir  ■'  i^'PQT       ^''i 
IT'I^PCT       ^6" 


1  ■"  '""J  ""^  t^ 

^-  r-Q  ~5  4 

=  '■  '-0  -■:'  ~ 

"•  ic  :g  ~ 


;  09 

:"ri.-i 


51-1 
il5 


f^'J       0    0  '     oO    '0    3 


- ; r-ENs :  crj  .^  :■  ■  .1 .  - . ,  7\  3 ,  i" 

CO    5    I--1.  3         ' 

' '  t ,  J !  -:.  0 

■r-  I     E^     j;    -■  :.  .,;  =  .   0 

;r|A^'/.   K>  -J.;  0     0       ;0      70      10 

DC    3    I=ivi,  J 

:  ^ .'  i  J  ( : ,  >. ; 

:0    7     J=i^,  : 
S!=A^i  ;<,  ^. 

-••■jv :,  j:'=3i 

7     ■nN-:M,.;=- 
1-0    TO    1  3 

;o  TO  z~ 
10  e :  =".j  ■;>.>;• 
:a  i:    ;f=-<.  = 

1  1    A.^'  ( ;<  .  •  ■-  ,  =i,.   :■:/.,  :  p  ;  ,  31 
tO    15    '=1,2 

1=^'    I         ::.  J        A  )      i"     -TQ     1  .-. 

3 1  ='»o  ',;,!*> 

00    1;     .=^,  :, 
1 2    A  .'  t  J .  J  :  =,^ ..  '  :  ,   .  1  -3  :  S..1,, ;  .•  u  .    . ' 

li   :crj7:f:i£  "■   ' 

Ic    CCNT.-iV^ 

"0  1  ■*  :  =.  1 ,  3 

CO     17    J= 1 ,   ■ 

17   )■■  t,  v>="j'  r.  ..;-3i 

■■JO    77     20 

S7  •...f^i'^zt-  i:3o-    ;- 
ic;c   PO^-iAT    .i,:_7;;;^   -,,^    ivi^,^:p<3t   c;CES   not 

tJ    Pt7'JPflJ 


■yt  r£-. 


■-^0 


DO    10    ''^l, '■  J   ."■.!"  0^7  "'     ''  ' 


164 


in  -t 

z'ii, 

242 
243 


247 

248 

2S0 


izb 


2.il 
i64 

c'-'C 

267 
S71 


;-'5 


^=•6 

oc-> 

5=5 

29'^ 

jiJO 

iSi 

■^Oj 

31 0 


C  I  J  )  =A  ,  v;, ,  J) 

■i-iA,  J I  =->  (  f,,  J  '  'A  I  ;■;,  J  •, 

CO    ''O    L  =  ^.  '3ivT^- 

Tr'NUNCDS    ,  L~        I  '^      "JO    ',0    20 

£o    40    LL=L.   in^jTM 

J-J+1 

A'  ;,  ..M=.i>(  I.  .J     -C.L:  tA(K,  LL) 
4C    C3NTI'J'0E 

-(  '  )  =-  :  '.'-'.'  L  '-^F'"/^.' 
3C    CCNTIMUf-' 
10    CCNTINI.e 
ICC    K=i<-1 

tF(K       EO       0^    QO    TC    2C0 

DC    50    J=2.  I3W7H 

IF'NU.NODS       LT       L"     C-O    TO    SO 

r  1  :■<  ;  =C  .'  ^  1  _^  1  ^  _   ,  r  ;  ,p  I  (_  ) 

50    CCNrilsiijE 
GO    TO    too 

200  ^C.MTINOE 
aO    C2NT:nlE 

201  -ETl.RN 
END 

3USf*CUT!N£    I3C'/E?i  \HX,  VXY:  Zl'VZ,  VOI  'JME ; 

,^',!:f''~^§^-^'-'^ ■---'■  ■^•'-'■-  '^•j(3.  61 .  -.-3;  31,  ■ 

2CNX(20),  'Ny<20^    CNTirO)  -^  » -- '  ■  t- -■  ■ --  J  • 

COMf-iQN    ::j,  -F,  3  I",  V  /    rr 
:3P0T=1 
iLPHA=0    57^35027 

/[=_ ALPHA 

■;g  to  10 

2  (  :=ALrHA 
'30    T^j    ;  C 

•JG    "0    ;■:■ 

4     Xl--.\\_PU^ 

;3  TO  10 
?     ''t=-ALrSA 

:iv.AL.--MA 

'iO  -u  1: 
=    .x;=Ai_F'--A 

'30  *0  10 
7     VT=ALrHA 

>e  Tc  i-: 

3  ,t  I^— -.L?-A 
10    XI==1     '^<' 

xri=:.  0-;: 
■-:?=■.  O'v: 
v::f=..  j-v; 
1 1  p  =  •   0  + ;  I 

CMi'I  '  :  >  =  ^r,--:-i.z;ri-*i  2.  o»>;  :^v;  +  7Xf  i    -■.',^    - 
ONX  I  (  2  J  =^  I'"»Z  ;  ^>  ■  2    0-*'<  J  --T -z  i-i ' 
DT-KI  '-  3  )  =v  "P-tZ  ;.'^,f  ■  2    i-jfi,  I  j.vi_7  r  _  ) 
CNX  [  <  A  ;  =/ rp.,  J  r  .<»  ,  .-T    .^>;(  I  _v  ;  *~  r.1.  f 

'-Tix I  (5'=v :.'■'■>: :'--*i 2  ''"ri-vi-z"  -i 

CN,"  ;  <i  !  =v  if^5^-  2--»  1  2  O'X  r-v  I  -Z  1-  ' 
L'N  X  I  <7  )  sr'iS'f;:?*  I'J  .-lej'i-',-  .Jt-T 
"NX.  :  tS.' ==•' ;,---ZI"->  •  2    .T^  (  J-,- r_  r '.  {      ;i.3 

Z fi V  r  ■•  -^ )  =-  - 1  r- » i  I  ;"■♦  / :  /  7  0 
z.^iX!  V :  l;)  =  ; .    :-vi>v;  >iz:M/4  o 

Z  N'  X  1  ;  1  ;  /  =-  '  ■  r  -»  ■'  T  ■■;■».  jf  r   -  T     -1 

CNXI  i.  2-.  =-ZN<:  •  iz  ;    '    ~   '■ 

zcj  <  r  <  1 3 ;  =  -  ■  1 ,  ""  -  z ;  ♦  z  i  ' » ■■' ;  r^  /  4  r. 

r.i-^xi  ( 14  ■  =-zMx: '  I  j^ 

r:^  X I  ( 1  = '  ■'  <  :  0  -  z :  -'  z  L .'  * V  i  p ,'  4  o 

..  f  J  X  i '  1  &  ■=  -  r,  i'  J .» .' '  ■.  •■■■ ' 

3-ixi  f  ;7'=--'i.-ifz;>^->\:   :■  .-■ 

rNxrii3:  =  ^*.   ■:-■•  :-M  T  ■.  g.r  rp  .  a  o 

ZNxr  f  i  =  .  =--v:p-.::^>Ai-  -2  -• 

ZNXl;20)-=-D-i\:  ,15) 

ZWVI  -r  •  ;  •=  x  '  .•" ^  Z  I. "-->■.  2    'T-f  ^  '  f  .< '  ■►Z  1  ■-  i 

;;?(■':■, 2 -  =  v  I,- *7  f-i»(;   •  » ,  r_,  ?  ^^t^  . 
:ny  I  ■  3  •  - ■;  r  =-*  z :  ■■•♦ ;  2  0*  v :  -n- 1  -  z  i  -  i 


■  '3 

/3 

y 
^ 

0 

,-', 

/  '^ 

0 

\J 

"8 

■J 

.", 

/'5, 

0 

■^  ^ 

-Q 

3 

n^jv  r  (  J.  1  =x  I  '1*Z  :.-i-; 


li«'/X- 


"  )  '5    0 

J ' .  s  ■: 

J  '  /6.  0 


-  ^  1  -  1    0  .'  .•■  3 


165 


31=. 
313 

|l5 

317 

31  = 
3S0 


3^:3 

330 
33 1 
332 
23-3 
-334 
335 
334 

339 
340 
3-it 
3^2 
3-i3 

B-j 

345 


"i3 
■350 


354 


353 
-^o 

360 
oil 

3ib3 
ji-l 
■365 

1*7 
36S 
36" 
3'0 
331 

373 
37  ** 
375- 


373 
37' 
330 
331 
362 
3G3 
'J?4 


335 
3£o 
337 
•?PS 
389 
390 
3-3; 
392 
393 


■J  )  /  i 

•"■ ;  ■'  p 


n  )  ..-  o 

0 

0).9. 

0 

0)/8. 

0 

CO/ 5. 

'- 

0)/0. 

'^ 

0)  .'i 

0 

o)/e. 

0 

0^'8. 

rt 

EN VI I  ij '■--."  I" :!,-'>•  a:  :.■'>' '-  <  ■ -1:-: 
DNYi  i6)=.<  ;?•>;.;-■»'.;  o»'  i-:'i-Z!+ ■ . 

DN'r'  \  <"  '}  ■-  ^  i  -  >Z  1-  ■> '  7    ■"•■♦■'■+  X'  1 1-  Z  !  -  1 

3riv  1  is.i  =,x  '.:-:»;;  =  »■. 2.  O'*'- 1  -  <i  *-z  :-i, 

DNV I  1'  ^  i  =-  (.  i    C- i'  I  *  X  i  *  ♦  "  I  '^  ' a    0 

c^i V I  u  c ;  =- ' :  t-'  ■» ." :  .'"•»  v  i  -'i  o 

Z'l'tl  ill"  =-DNY' I  ■■  '=  i 

Df-iV  I  1  1  3)  =-v  !     ."i-Z  :*2I  )»  i  r^.a    0 
DNVI  (  14;=-(  1    0-Z:*ZI  >-»UP/4.  0 

DNv  I  <  1 5  •«  = '  ; ,  e -: :  > :  1 1  * ^'  r p  ,■  4  o 

CNt'I  (  Ib)=-CNv!  V  13) 

dmy:  (i~'---'..  ;■- V 1 -t^i  i  +  zip/'i,  o 
DNYi  ( 13"=-';:=  ^zip-y: /2  0 

Dfivi  ,'  i  =  -=-z,Nv:  :  \7; 

3.NiV  I  (  20  J  =-  '  :  'n* Z  I "-» Y I  '2    C 

DNZI  <  1  ;=Klri*^I"!»<r    OfcZ:-<I  i-r  lJ-1. 

ONZ  I  r5  .'=x  ;P>v  !ri»',  2   0-tZi-X  I-^YI  +  i 

DNzr f3!  =  if :p-»f iP-n's  c*z:-<r-Yi-i. 

DNZI  C4)=.«,  I^-»v  !p^<^    0>;'T  -XI-YI*! 

DMZi(  5)  =  «  it<>^-r-'^<  J.  o»zi-<:  -''I  -1 
DNZI  («:j  =  x  ;p  >,  ir*  :^.  :->zi-^xi-yi-i 
ZMz:  (~)=:f  ;p-»vip>,  2  0*7  r*'<r*vi-i 

3MZI  (31  =X  .".'^■»'--  !r»'  3    0-i-ZI  -<i*YI-l 
DN  Z I  ■''>=-■  1    .-,-<;*<;,*  r  I M  /■  4   Ti 
DNZI  '  10;=-'.  :     Z-r  !-*<:  ;^XIP,.  4    C 

DN z  I ;  1 1  ■■  =-  ■  1 .  ■:•-  X I  -t  X  i  ■;  *  / 1  p  /  4  0 
CNzi  ( 12)=-.  1,  :--^':  +  -  i  !*x:r!/4  o 

D^JZI  (  I3>=-XI'"^Y*^"■*ZI- 2    D 

DMZ I ( 1 4 ) =-x iPtY :m»z I / 2  0 

DNZ  I  ■;  1  5  >  =-''  rF-t'TP-fZ  ;  -'2.  o 
DNZ  I  I  16!  =-X  r.^ov  iP*Zi,  2    0 

CNz:  ( i7)=-CNz:   ?:■ 

3NZ!  (  lS!=-t;NZ;  (10) 

DNzr  ^-  i  =  )=-rNz:  ■  1 1  > 
DNZI  ;2o;=-c'iZ:   ;2) 

TA.J==0    0 

Du  Ir  1=1.  iO 
1  5    TAJ=  T A  J-^DN  >•  I  ;  r  /  -*e  X  '  I  ! 

AJ( 1, I '=TaJ 

'■A.J=0    0 

□0  16  1=1, 20 
1  6  TA.J  =  TAJ-rr;N<I  ' 

A-.'  1.  2,'  =  TAJ 

TA.J=0.  0 

DO  17  :  =  1,  20 
17  TA-.^  =  TA^>*-D''J '  I  • 

Aj(  1,  3;  =Ti,j 

-AJ=0  0 

DO  13  1=1. 20 
1  Z     7AJ=TA.^-r[;NV  II':  -iFi  (  I  i 

A.^(  2,  1  ;  =TAJ 

TA„  =  0.  0 

D'3  19  :  =  1 ,  20 
i^  ~A^=7A..V0MYI  ■  I  ■)-*£•-{  I  1 

AJ(  2.  I ;  =^"-0 

'A,^=0  0 

DC  20  :=' , 20 

20  TAv'=TA...-CNx;  >  I  :  +ez  ■;  I  ) 
Aj'2. 3  =7Aj 
"A.J=0  0 
DO  21  1=1.  20 

21  ~a.^  =  ta;»c.''JZ  r  '  '  •  -<i=^x  (  I  '' 
A.^3.  1  ;=Ta„ 
7A,J=0  0 
DO  22  :  =  1  20 

22  T.^j!--JA.^-rC-IZ  I  '  I  '  *£Y  •  I  ■■ 

A.J.  3,  2.='A,., 

'"Aj=0  0 

DC  23  1  =  1 .  -:o 

23  T«.v=TA.J-D■■^iZ  ■  i  1  ■  »f  Z  <  I  ' 
AJf  3.  3)=TA-.- 
OJ=»A.;  f  t  ■  '  >  "  Aj  ,  .2,  r  .  i-Aj  '1,3;  -AJ  (  !  ,  1  >  i-4J 

i*Aj(  1.2;  »'■•-'• ;;  :>  '-'a  j  "i,  i  •  -aj'  -. ,  2 )->'..  <  2 

2 "-AJ'  1.  3  t  "■■■."  2.  1  .  *A.  '  J,  n  )  -A.J'  1 ,  3  >  >A  J  (  2 
Ci^AX=A35;D.  ■ 
ZALL  IMV.E.P^  ;  -.J.  r.  I£i 
DO  30  I-  [  .  D'l' 

DNx  ;  I ;  =-- ^  1  •  '  ''C'..  < !  ( ;  1  +  r :  1 .  2 ;  ^Dr-jv  • .  i  .  ■■ 

r-.tjv  :  I  )  i^T  .;  2  1  •  *ZMX  J  ■.  I  i  -^i  ■■  2.  2  J  »0,">'>  :  I  >  - 
30  CNZ  .  r  ,  =T  ;  :,  1  ■  »,,NX  I  (  :  )  <-T'  3.  2  '  ■>r.N-  I  '  ;  >  • 

DO  j5  I = ; ■ 20 

CO  35  j=l. 20 
35  EA'  1,  J  ■  -EA  I .  j;  *  ■■  (y  X»DNx  C  I  ;  'ONX  .  J  i  *v'a  X-DNV 


*£■'  ^  I  ) 


:  > 


■»■  A  ..;  (  3  . 


,  1  '  »A.. 
■  ?  )  ■<  A 


Ti  , 


:/■'':.  N  Z  I 

!  .■  »i;-rjZ  I 
i)fONZI 


•  DN 
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3  =  4 

-■  J 

-00 
10  1 
405 
an3 
4C4 


V'CLUME=Vi3L'.'r^E*n.jAX 


IF( ISPGT 

ir 

2 ; 

C-G 

'19 

2 

if; :srGT 

-     1.   J 

■^  \ 

-;o 

3 

^  '^   i      ♦  '-  ~  f-.  T 

-  ■    \  A  c  r  1-  ■ 

r.  J 

.;g 

!  U 

IF'  ISF'^T 

-  "j 

^  ■■ 

C-C 

Tn 

:^ 

IF( ISPGT 

io: 

3  ) 

50 

"^,-n 

'o 

I-CI3PC-!- 

^0. 

7) 

■5G 

TG 

IF' ISPG7 

E>3, 

3 '. 

50 

TO 

3 

^ETURfJ 

£ND 

"YPt, 


406 

40  ■' 
4  OS 
40'^ 

4:o 
4: : 

413 

4;  4 

4  1o 
4t- 

4 ;  5 


4^3 
430 

4  32 
-133 
4->a 
435 

4''.*5 


4^0 

441 

44S 

443 

.144 

44  = 

446 

447 

44.3 

449 

450 

451 

452 

453 

454 

455 

4^6 

457 

453 

459 

440 

4^1 

462 

463 

4.!,5 
4t4 
447 
4=3 
4«,f 
470 
471 
472 
473 


■rUERGUTINt-:  L3CA'  (riUNODS-  NUMELS.  ^' ,   OG.  '.  ^,  ' 

D I  MENS  I GN  .'  <  ;C0  ':.^".  ZOO  ).Z\  SCO  ■  ,  Th  ( 3  .■  .  -l(  =0-  20  ^  ,  5NV  (  200  i 
1.  i  I  i-FEi  3<.'  ■  .  .3:Ai^f  50  ;■,  PN(200)  .  THNt200> 
ThA=0.  3235998 
03F=GQ/3.  0 
■Sl2'^=— 3Q/;2.  0 
ST-SINi  TPA  ,. 
3T2=SIN(THA/2  0) 

:'=cgs.-tha; 

3T2=003' THA,-2  3,' 

J.J=0  0 

ffEAD  '.r.  lOlO.'  R,  TH<1)  ru(2),T4(3) 


0.  CN'/,  r.i, 
. a (200; 


;r.  TMi  i ) .  TH'; 


5  -E-C  .-5,  :oi 
=  DO  10  :=i, 3 

KN'  J.J;=SP 
THiN  I  JJ  ■'  =0  0 

r  (  :-j )  =0  0 

f  (  J J  i  =c 

z ( vj)=TW( r  i 

10  IGi-jriN'.jE 

=N.' Jj')=R9 
rHNi  Jvl=Ti-iA/2.  0 
X  (  J  J)  ■=R»3T2 
r'.'J^'-,  =p-»CT2 
Z'' JJ;=TH^  1  ) 
■Jj=J..-'l 

THN-'-JJ  i=ThA/ ~  0 
X(  JJ-=R-»ST2 
V(.JJ)=R*CT2 
Z  >  Jij  ,1  =  T-li'  3  ' 
DO  20  t=l, 3 

FN,-  JJ-,  ==R 

<  '  Jj;=^*'^T 
r  (J.j)=o»CT 
J(.vJ)=THf  i  ) 
20  -.GNTIfJuE 

If  (v'-j  .  oe     ■■^ijnod3>    ^o  "0  100 

PtADtS.   1020>RR,  T'r>(  1  )     TH(C^) 

DC    3-0    1  =  1    2 
-v=.- J-1 

T-ri(  .J,;  i  =0     ) 

A    '     J-J    I   =0        0 

30  :'■  jji=TH( : , 

DO   40    r  =  i , 2 

JJ  =  OvJ+l 

Rn; jj)=^r 

■^hM,'  JJ)  =THA 
X  ,  JJi=R-t3T 
/'.  J.J)=P<>-3  7 
40    7  'v>.';=rH<  I  ; 
■50    TO    5 
iOO    NUl ^MUNCD^-e 

DC    101     I*.     -j';i 

101  a,\v(  [  ■=-!    ; 

."JIJl^NUNCDS-" 

:0    102     I-fjiJl    Ml.'NGDS 

102  RNV(  I  1=0    0 
J  '  1  )  -"'J  1  .tM 
li  (  2  ■- ■'•■J3P 
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'l^:?  G(4)=i33P 

«""'-'  G<3>=0i;(^ 

"50  DO  i05  :=9.  »it'NGD5 

^Sl  105  3( : 1=0    0 

-32  i  10  Nil,  1  i'^l 

^82  Ml  1,  2) =3 

"3-1  Ni  i,  31=^0 

435  Mfl,'i;  =  12 

iSi  N  (  1 ,  5  ;  ='  1 

'-c.a  N<  ! ,  7  ;  =  1  = 

^39  ^(',  1,  3)  =  13 

4^0  M( 1. 9)=5 

-J'^l  N<  1,  10)=12 

J^i  N( 1,  i  1  ;  =  i7 

4=?"  ^J^  1,  12)  =  !G 

4-4  N( 1. 12  1=2 

-^5  .■\J(  1,  14. =7 

-'^^i  N(  1,  lr/  =  I  = 

4  =  "'  ?J(  1,  16)  =  14 

4-3  N(  1,  17-.  =4 

N  (  1 ,  1  3  )  =  !  1 

N  <  1 ,  1  =  )  =  1  6 

N( I,  20 > =9 

502  I"'YPE'1/  =  1 

-'0."  LEAK'.  !;=1 

292  -■"  1-^'^  :E=2.  NUr-^ELS 

305  DO  120  J=l, 20 

50o  12C  NJ  '  ic.  J>=r-J(  lE-l,  J)'12 

^07  ITVPE':  !E)  =  l 

505  L£AK(IE;=1 

509  lao  CCNTIi'JUE 

510  lOiO  FnsMA7(4Fio  5) 

511  1:20  =^0R^'Ar;3Fi0  5) 
rii  SETURfJ 

513  END 

Z--^  S'-'SPOUTINE  ■i,.''jSi"':tNUNODS,  .?T,  X,  Y,  DEP'^K-  F,  N,  PN,  -HN,  ,■-" 

i  ;  -  pi "^'-'^  i  jl'J  N  .-  30,  20  '  .  F  '  200  ; ,  RN <  20C  :■ .  THN  (  200  )  .  X  '  200  ) ,  /  '  PCO  ) 

- 1  i  ^^i^J  '.it.  1  ) 

51-  PT=' r  ■;  .: '^le?"'-!  i  ■»-9eri2  0 

518  PS  =  CEP'"H-»  10047    r\ 

-  t  9  ■-'  1  =N  (  M-  1 ,   1  ■■ 

520  P  1  =  ' F(  .^1  ,' *DEP"rH  ! -*ogo2    0 

521  3LCP=i -=7-^;  >  ,  (RT*(  1    0-PI\i(Jl))) 
5-2  =;T=PT*0.  5-»iPS— =T)/3L0P 

5^-  CO    10    i=9,  •■JUNGuS 

524  TH=THM(  I  .' 

525  <i  I  ■=SN(  I  ;*r;T»SI"J'TW) 
^-=  V  '  I  l=RrJ(  I  '*PT*CC3'  TM) 
?2~  10  CONTINUE 

523  PET'JPfJ 

52°  END 

-^^'  SUePCUTINc  A.^i.Sf'N,  M,  NUMELS.  Z.  7FPTH.  H  ) 

-31  D I rENS ■ GN  N •  50 .  20 ) •  2 i 200  ^  F ( 20C  ' 

rj2.  CO  40  I=ri.  NUrELS 

5;3  DO  10  I !=3. 4 

534  .■=N(  r,  r  I . 

|35  :Z=41.  0031t.>LEPTM-40.  0-»  i  F  (  J  WDEPTH  ) 

53o  ;F'ZZ       LT      O    0'-     ZZ  =  Z  ■  J '■•  ' "    " 

°3Z'  .-     IF'::       OE.      OEPTM)'"zZ  =  iOEPTH-Z(  Jj  ) /2.  O+Z'.J) 

3jS  10  :'j'=zz 

539  :a  15  :i  =  io.  12 

540  ^  =  N(  I ,   III 

|-*1  ZZ=4I.  C0814»[;E=7m_4o.  04(F' J)+DEPTH) 

5-13  IF(ZZ      LT      0    0;     :Z  =  Z(J>  -2.  0 

i'\-i  _     :i-<ZZ       i,t       DE=Tm1     ZZ  =  (DEP''-H-Z(  J)  )/2.  0+Z(  J) 

-«>4  1-    i(J)  =  ZZ 

343  DO   20    r:  =  15,  1  = 

544  ■j'T-j,-  :.   I  ;  ! 

547  .^■=N'  I,  I  r■-l2.• 
5■ie  20    Z  '  J  ■  -■• '  CEP  TH  .•::..  '.J  ■)  /r   :-i 
54=-  40    CONTI'-lcE                                "■   ' 
■550  'iEa'PN 
5  5 1  END 

SEnTRv 

JOB    FAPArE~E,9i; 
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